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Figure 3.46 Cell perimeter; (a) the normalised mean, i.e. (sij)/p(s11), and (b) the standard 


deviation, o(s;;)/o(811). 


: 


1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 
1.0 1.0000 0.9376 0.9227 0.8878 0.8918 0.8717 0.8324 0.8271 0.7906 0.7835 
0.9 0.9437 0.9333 0.9039 0.8577 0.8241 0.8046 0.7922 0.7906 0.7556 0.7472 
0.8 0.8996 0.9127 0.8580 0.8266 0.8131 0.7558 0.7752 0.7560 0.7580 0.7441 
0.7 0.8840 0.8766 0.8334 0.7928 0.7644 0.7313 0.7054 0.7042 0.6937 0.6788 
0.6 0.8780 0.8339 0.8473 0.7684 0.7392 0.7261 0.6757 0.6806 0.6575 0.6363 
ry 0.5 0.8509 0.7948 0.7592 0.7228 0.7301 0.7017 0.6550 0.6364 0.6166 0.6036 
0.4 0.8019 0.7457 0.7649 0.7284 0.6798 0.6613 0.6011 0.6036 0.5858 0.5847 
0.38 0.7852 0.7609 0.7129 0.6883 0.6618 0.6408 0.6161 0.5767 0.5755 0.5662 
0.2 0.7860 0.7273 0.7305 0.7008 0.6431 0.6189 0.5638 0.5476 0.5334 0.5264 
0.14 0.7771 0.7198 0.6895 0.6911 0.6296 0.6090 0.5538 0.5440 0.5124 0.5003 


Table 3.22 Numerical means of the cell perimeter of compressed Voronoi. 


Figure 3.47 shows the plot of the change in the perimeter of faces when the Voronoi structure 
is compressed. The normalised mean and standard deviation are shown in the form of contours. 
The values before normalisation are p(u(zi;)) = 0.3766 and o(pu(z;)) = 0.0532. 


Both the surface and the perimeter are embedded in three dimensions, the former as 2-d facets 
while the latter as 1-d facets. But, comparing Figure 3.46 (a) with Figure 3.46 (a), the trend of 
changes in the mean is smoother for surfaces than for perimeters. The difference between the thrend 
of the standard deviation is even more pronounced, as can be readily seen by comparing Figure 
3.46 (b) with Figure 3.46 (b). This may imply the reduction in the correlation between different 
statistical properties as the structural differences increase, as surface is directly related to the convex 
hull of the volume while perimeter is a convex hull embedded in another convex hull. Also, this 
means that convex hull as a mapping or function is not smooth, or it could mean that it is a smooth 
function only up to the order one, that is to say, when we apply it only once not twice or more. 
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Figure 3.47 Face perimeter; (a) its normalised mean, i.e. (u(zj))/u(u(z11)), and (b) the 
standard deviation, o(u(si))/o(u(s11))- 


Ve 
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 

1.0 1.0000 0.9640 0.9460 0.9103 0.8907 0.8659 0.8523 0.8339 0.8155 0.8056 
0.9 0.9620 0.9384 0.9040 0.8681 0.8371 0.8186 0.8124 0.7927 0.7717 0.7606 
0.8 0.9262 0.9137 0.8738 0.8409 0.8129 0.7823 0.7755 0.7608 0.7489 0.7372 
0.7 0.8956 0.8722 0.8355 0.8058 0.7801 0.7378 0.7246 0.7276 0.7053 0.6874 
0.6 0.8745 0.8404 0.8214 0.7807 0.7527 0.7252 0.6936 0.6849 0.6661 0.6504 

Ty 0.5 0.8485 0.8093 0.7808 0.7514 0.7255 0.6977 0.6734 0.6479 0.6368 0.6283 
0.4 0.8186 0.7780 0.7614 0.7284 0.6951 0.6659 0.6335 0.6139 0.6065 0.5987 
0.3 0.8125 0.7771 0.7347 0.7075 0.6670 0.6490 0.6210 0.5916 0.5772 0.5731 
0.2 0.7958 0.7459 0.7225 0.6946 0.6506 0.6254 0.5870 0.5685 0.5528 0.5466 
0.1 0.7938 0.7451 0.7154 0.6897 0.6433 0.6181 0.5815 0.5594 0.5350 0.5196 

Table 3.23 Numerical mean values of the face perimeter in a compressed Voronoi. 
o (u(z4j)) ee 
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 

1.0 1.0000 0.9795 0.9413 0.9524 0.7858 0.9426 0.9614 0.9060 1.0364 1.0311 
0.9 0.9221 0.8227 0.8754 0.9145 0.8047 0.8613 0.8197 0.8653 1.0369 0.9671 
0.8 0.8363 0.8407 0.7857 0.8670 0.7517 0.8090 0.8055 0.8980 0.9424 0.8702 
0.7 0.8939 0.8322 0.7718 0.7612 0.7903 0.7904 0.8453 0.8296 0.9233 0.9350 
0.6 0.8724 0.8413 0.8031 0.8155 0.7560 0.7662 0.7972 0.6784 0.8740 0.8279 

Ty 0.5 0.8689 0.8957 0.8640 0.7293 0.6592 0.7756 0.8340 0.7986 0.8512 0.8750 
0.4 0.9538 0.8608 0.8199 0.7414 0.7733 0.7324 0.7916 0.7954 0.7916 0.8463 
0.3 0.9359 0.8457 0.9174 0.7016 0.7422 0.7711 0.7139 0.7697 0.7137 0.7206 
0.2 0.8985 0.8447 0.7545 0.8238 0.8155 0.7855 0.7918 0.7724 0.8133 0.8138 
0.1 1.0322 0.8763 0.8181 0.7953 0.6991 0.7159 0.7692 0.7085 0.7546 0.8212 


Table 3.24 Numerical standard deviations 


of the face perimeter in Voronoi network under 


various degrees of compression. 


The simulation above is then repeated again for a second time, where now p4(s11) = 2.1206, 
o(811) = 0.8154, w(u(z11)) = 0.3840 and o(p(211)) = 0.0518. Figure 3.48 shows the results for both 
the cell perimeter s and the face perimeter z. Also, compare these pictures with Table 3.25. 
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Figure 3.48 The cell perimeter and the face perimeter; (a) (si)/p(si1), (b) o(sij)/o(s11), 
(c) m(e (zig) /MCu(211)) and (d) o(m(zig))/o(u(r1)) 


Ty 


0.9 0.8 0.7 
0.9781 0.9164 0.9188 
0.9292 0.9004 0.8959 
0.9185 0.8512 0.8339 
0.8868 0.8596 0.8327 
0.8621 0.8241 0.7817 
0.8281 0.7785 0.7604 
0.7661 0.7653 0.7156 
0.7999 0.7731 0.7409 
0.7491 0.7368 0.6950 
0.7332 0.7161 0.6898 


Te 
0.6 
0.8766 
0.8330 
0.8090 
0.7596 
0.7425 
0.7350 
0.6856 
0.6744 
0.6551 
0.6345 


(a) 


0.5 
0.8417 
0.8144 
0.7808 
0.7573 
0.7278 
0.6821 
0.6833 
0.6747 
0.6370 
0.6265 


0.4 0.3 0.2 0.1 
0.8126 0.8282 0.7834 0.7633 
0.7851 0.7555 0.7511 0.7323 
0.7475 0.7256 0.7120 0.7015 
0.7262 0.6988 0.6964 0.6763 
0.6810 0.6723 0.6474 0.6397 
0.6752 0.6475 0.6352 0.6116 
0.6423 0.6052 0.5874 0.5698 
0.6093 0.6080 0.5627 0.5641 
0.6231 0.5672 0.5677 0.5361 
0.5977 0.5548 0.5552 0.5290 
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re 
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 

1.0 1.0000 0.9773 0.7894 0.9984 1.0195 0.8271 0.8252 0.9487 0.7626 0.7128 
0.9 1.1772 0.9638 0.8714 0.9202 0.8793 0.8753 0.8798 0.7462 0.7780 0.6670 
0.8 1.0766 0.9053 0.7051 0.8686 0.7943 0.9349 0.5831 0.6456 0.6093 0.6356 
0.7 1.0506 1.0094 1.0137 0.8884 0.7487 0.7961 0.6260 0.7008 0.7110 0.6383 
0.6 0.8534 0.8277 0.8383 0.7692 0.7554 0.7554 0.6306 0.6991 0.6439 0.5835 

ry 0.5 0.9281 0.8727 0.7584 0.8685 0.8945 0.7790 0.8540 0.8201 0.7273 0.6109 
0.4 1.0266 0.7761 0.8214 0.7686 0.7893 0.7270 0.76385 0.6527 0.6022 0.5956 
0.3 1.0474 0.9116 0.7985 0.7823 0.6053 0.8193 0.6155 0.7399 0.5585 0.5668 
0.2 0.8807 0.8048 0.7829 0.7438 0.6276 0.6633 0.79386 0.6742 0.6638 0.6000 
0.1 0.8996 0.8820 0.8161 0.8172 0.7064 0.6345 0.7206 0.6647 0.7240 0.6899 


(b) 
BHC 265 )) ee 


1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 
1.0 1.0000 0.9740 0.9327 0.9028 0.8710 0.8447 0.8119 0.8070 0.7864 0.7749 
0.9 0.9817 0.9287 0.9023 0.8735 0.8320 0.8126 0.7782 0.7608 0.7442 0.7353 
0.8 0.9482 0.9068 0.8681 0.8373 0.8039 0.7809 0.7508 0.7248 0.7079 0.7004 
0.7 0.9246 0.8814 0.8465 0.8174 0.7699 0.7500 0.7208 0.6933 0.6843 0.6670 
0.6 0.8846 0.8544 0.8126 0.7773 0.7453 0.7189 0.6878 0.6667 0.6467 0.6424 
Ty 0.5 0.8510 0.8268 0.7861 0.7496 0.7240 0.6879 0.6644 0.6370 0.6243 0.6118 
0.4 0.8542 0.7924 0.7657 0.7184 0.6930 0.6720 0.6354 0.6071 0.5865 0.5697 
0.3 0.8378 0.7876 0.7612 0.7213 0.6755 0.6554 0.6135 0.5900 0.5643 0.5549 
0.2 0.8126 0.7632 0.7354 0.6962 0.6575 0.6289 0.6062 0.5653 0.5530 0.5329 
0.1 0.8088 0.7494 0.7170 0.6906 0.6449 0.6155 0.5873 0.5531 0.5413 0.5172 


(c) 
r 


1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 
1.0 1.0000 0.9998 0.8469 0.9409 0.9015 0.8691 0.9077 0.8331 0.8906 0.8925 
0.9 1.0290 0.9684 0.7922 0.8566 0.7914 0.8284 0.8495 0.7599 0.8379 0.8420 
0.8 1.0285 0.9118 0.8258 0.8582 0.7783 0.8267 0.7019 0.7007 0.7504 0.7876 
0.7 0.9241 0.9131 0.9166 0.8512 0.7612 0.7458 0.7123 0.7265 0.7861 0.7844 
0.6 0.9830 0.8897 0.8278 0.6931 0.7288 0.7882 0.7028 0.6968 0.7300 0.7506 
Ty 0.5 0.9891 0.9817 0.7976 0.7720 0.8001 0.7218 0.7600 0.7386 0.7329 0.7131 
0.4 1.0478 0.8275 0.8022 0.7546 0.7731 0.6930 0.6875 0.6249 0.6975 0.7238 
0.3 0.9820 0.8738 0.8755 0.7251 0.6961 0.6874 0.5638 0.6426 0.5881 0.5863 
0.2 0.9251 0.7740 0.7667 0.7174 0.6943 0.6578 0.6800 0.5811 0.6519 0.6929 
0.1 0.9230 0.8312 0.7762 0.7384 0.7160 0.7186 0.7152 0.6357 0.6987 0.7378 


(d) 
Table 3.25 Numerical values of cell and face perimeters. (a) w(sij)/p(s11), (b) o(si3)/o(s11), 
(c) m(e(zij))/MCu(211)) and (d) o(m(sig))/o(u(s11))- 

Codes for finding volume, surface area, cell- and face perimeters can be found in § A.28. There 
only some of the printing commands have been left out. 

The procedure in general follows Algorithm 3.3. The volume is the summation of all tetrahedral 
volumes obtained from the triangulation. The surface area is summed over all triangular faces 
of a convex hull. Plane parameters are then calculated for all of these faces, viz. a = |1, y;, zl, 
b = |x;, 1, z;|, c= |x, yz, 1| and d = |x;, y;, z;|, in order to group them together into polyhedral faces. 


Algorithm 3.3 Volume, area and perimeter algorithms. 


(v, Ue) + Voronoi tessellation; 
exclude boundary vertices; 
exclude boundary cells; 
for every cell do 
find its Delaunay triangulation; 
V & YO | 25, ys, 22,1] /65 
find its convex hull; 
At ¥i(s(s — a)(s — b)(s — 0))/?; 
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for all facets of the convex hull do 
(a,b,c, d) + plane parameters; 
group hull facets into polyhedral faces; 
for all hull facets in every face do 
{vy} < vertices of all coplanar facets; 
find convex hull of uy; 
B < sum area of these hull segments; 
count number of occurrences of their edges; 
{er} < edges counted only once; 
{ec} + ef3 
2 D(D(Aai)?)!?; 
endfor 
endfor 
8 & (Dai)/2; 


endfor o 


The program varea.m was written long before vareac.m. It is listed after the latter in § A.28, 
although most of the variable names have been changed. The old names were long, for example 
CubeNormalVolumePerInnerCell1 which is now changed into cbnafn. The original program has not 
been published (Tiyapan, 2001, KNT8()). 

Some of the past results are shown in Figure 3.49. 


(a4) (b4) (e4) 
Figure 3.49 Past results of area, volume and perimeter. (a), (b), (c), (d) and (e) are respec- 
tively expected value, variance, variance of mean-normalised values, expected cube-normalised 
value and variance of the cube-normalised values, whereas (1) are values for surface area, (2) 
the volume of cell, (3) the area of face and (4) the perimeter of face. 


In addition, there is Table 3.26 lists the numerical statistics of the area and perimeter of face. 
Both this table and Figure 3.49 are from an unpublished work (Tiyapan, 2001, KNnT8(i), ibid.). 
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Ty 


1.0 
0.0322 
0.0301 
0.0280 
0.0261 
0.0242 
0.0224 
0.0206 
0.0191 
0.0177 
0.0166 


1.0 
0.9 
0.8 
0.7 
0.6 
0.5 
0.4 
0.3 
0.2 
0.1 


5(y] 


1.0 
0.6952 
0.6733 
0.6526 
0.6330 
0.6147 
0.5981 
0.5832 
0.5706 
0.5605 
0.5537 


0.9 
0.0300 
0.0280 
0.0260 
0.0241 
0.0223 
0.0205 
0.0189 
0.0174 
0.0161 
0.0150 


0.8 
0.0278 
0.0259 
0.0240 
0.0222 
0.0204 
0.0188 
0.0172 
0.0157 
0.0144 
0.0134 


0.9 
0.6728 
0.6504 
0.6291 
0.6090 
0.5903 
0.5731 
0.5578 
0.5447 
0.5342 
0.5270 


0.8 
0.6514 
0.6285 
0.6067 
0.5861 
0.5668 
0.5491 
0.5332 
0.5196 
0.5087 
0.5012 


0.7 
0.0257 
0.0239 
0.0221 
0.0203 
0.0186 
0.0170 
0.0155 
0.0140 
0.0128 
0.0118 


0.7 
0.6313 
0.6079 
0.5855 
0.5643 
0.5444 
0.5261 
0.5097 
0.4956 
0.4842 
0.4763 


Ve 
0.6 0.5 0.4 0.3 0.2 
0.0237 0.0218 0.0200 0.0183 0.0169 
0.0219 0.0201 0.0183 0.0167 0.0153 
0.0202 0.0184 0.0167 0.0151 0.0138 
0.0185 0.0168 0.0151 0.0136 0.0122 
0.0169 0.0152 0.0136 0.0121 0.0107 
0.0153 =0.0136 0.0120 0.0106 0.0092 
0.01388 0.0122 0.0106 0.0091 0.0078 
0.0124 0.0108 0.0092 0.0077 0.0064 
0.0112 0.0096 0.0080 0.0065 0.0051 
0.0102 0.0086 0.0070 0.0054 0.0039 
(a) 
le 
0.6 0.5 0.4 0.3 0.2 
0.6126 0.5956 0.5806 0.5679 0.5581 
0.5886 0.5711 0.5555 0.5423 0.5320 
0.5657 0.5475 0.53813 0.5175 0.5067 
0.5438 0.5251 0.5082 0.4938 0.4824 
0.5234 0.50389 0.4865 0.4714 0.4593 
0.5045 0.4844 0.4662 0.4504 0.4377 
0.4875 0.4667 0.4478 0.4313 0.4178 
0.4727 0.4514 0.4318 0.4145 0.4002 
0.4608 0.4389 0.4186 0.4006 0.3856 
0.4525 0.4301 0.4094 0.38908 0.3750 
(a) 


Table 3.26 (a) The area of face and (b) the perimeter of face. 


These much earlier results show in addition the face area statistics. Another set of codes to 
do the same job has been developed in the style presently used. It used to be ghull then, but now 
voronoin on Matlab is used instead. This is the reason why there are many programs listed in 
the appendix, from page 213 to page 304. It is a good practice to use more than one program for 
doing the same job, in order to cross check between the programs. As a program develops, the style 
and approach of the programs he writes also change. This is why it is difficult to leave out any 
certain piece of codes. What you see here are in fact merely snapshots that have survived from the 
continually changing virtual working ground. 

The results for the face area of a VT under compressions are made up of u(u(Gi;)), o(w(Bi;)), 
p(a(Giz)) and o(o(G;;)). These are normalised respectively by p(pu(Gi1)) = 0.0065, o(u(fi1)) = 
0.0018, w(o(611)) = 0.0065 and o(o(611)) = 0.0029, where 6;; is the face area under a compression 
such that r, = 7 and ry = j. These resusts are shown together as Figure 3.50. 


{ 


(u(face area)) 


W(o(face area)) 


0.1 
0.0158 
0.0142 
0.0127 
0.0112 
0.0097 
0.0082 
0.0067 
0.0052 
0.0038 
0.0025 


0.1 
0.5516 
0.5251 
0.4995 
0.4747 
0.4511 
0.4289 
0.4084 
0.3900 
0.3746 
0.3633 
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Figure 3.50 Statistics of face areas under compression; (a) the normalised mean of mean 


be(w(Ba3))/e(u(P11)), (b) standard deviation of mean o(u(fi;))/o(u(P11)), (c) mean of stan- 
dard deviation p(o(Bij))/m(o(B11)) and (d) o(a(Bi;))/o(o(611)) standard deviation of standard 


deviation. 
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Table 3.27 has the normalised numerical values used for plotting these contours. 


1.0 
1.0000 
0.9052 
0.8316 
0.7487 
0.7221 
0.6392 
0.5825 
0.5183 
0.4633 
0.4478 


1.0 
1.0000 
1.1387 
0.8035 
0.7238 
0.7988 
0.6963 
0.6974 
0.7246 
0.6779 
0.6377 


0.9 
0.9427 
0.8674 
0.7443 
0.7311 
0.6408 
0.5815 
0.5418 
0.4655 
0.4208 
0.3951 


0.9 
1.0757 
1.0314 
0.7713 
0.7318 
0.6993 
0.6241 
0.7181 
0.6884 
0.6656 
0.5708 


0.8 
0.8435 
0.7970 
0.6698 
0.6343 
0.5886 
0.5073 
0.4803 
0.4324 
0.3781 
0.3437 


0.8 
1.0285 
0.9593 
0.7159 
0.6668 
0.6248 
0.5705 
0.5614 
0.6228 
0.5349 
0.4823 


0.7 
0.7402 
0.7112 
0.6590 
0.5752 
0.5331 
0.4759 
0.4062 
0.3713 
0.3317 
0.3059 


0.7 
0.8041 
0.8929 
0.7835 
0.5928 
0.6091 
0.5704 
0.4507 
0.5051 
0.4740 
0.4569 


le 
0.6 
0.6837 
0.5991 
0.5745 
0.5438 
0.4699 
0.4408 
0.3724 
0.3183 
0.3141 
0.2757 


(a) 


le 
0.6 
0.7212 
0.5949 
0.6792 
0.5840 
0.5329 
0.5172 
0.4437 
0.4242 
0.4543 
0.3863 


(a) 


0.5 
0.6645 
0.5730 
0.5346 
0.4958 
0.4371 
0.3837 
0.3489 
0.2948 
0.2619 
0.2296 


0.5 
0.6669 
0.6301 
0.5606 
0.6665 
0.4589 
0.3629 
0.3370 
0.3712 
0.3349 
0.3285 


0.4 
0.5893 
0.5306 
0.4893 
0.4651 
0.3904 
0.3311 
0.2880 
0.2527 
0.2139 
0.1870 


0.4 
0.6257 
0.5234 
0.5531 
0.5509 
0.4261 
0.3178 
0.2802 
0.3013 
0.2668 
0.2744 


0.3 0.2 
0.5488 0.5002 
0.4914 0.4389 
0.4592 0.3924 
0.4474 0.3593 
0.3614 0.3302 
0.3157 0.2762 
0.2522 0.2293 
0.2132 0.1886 
0.1702 0.1422 
0.1410 0.1030 


0.3 0.2 
0.6937 0.5409 
0.6477 0.4421 
0.5506 0.3955 
0.6229 0.3910 
0.4609 0.3835 
0.4090 0.3335 
0.2536 0.2256 
0.2415 0.2232 
0.2204 0.1519 
0.1997 0.1332 


0.1 
0.4608 
0.4155 
0.3679 
0.3512 
0.2965 
0.2468 
0.2007 
0.1624 
0.1151 
0.0714 


0.1 
0.5482 
0.4561 
0.4186 
0.3941 
0.3318 
0.2881 
0.2190 
0.1792 
0.1118 
0.0677 
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. 
1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 
1.0 1.0000 0.9084 0.8333 0.7368 0.6714 0.6736 0.6061 0.5907 0.5819 0.5499 
0.9 0.8728 0.8531 0.7962 0.6763 0.5947 0.5654 0.5499 0.5301 0.5263 0.4965 
0.8 0.8034 0.7278 0.6479 0.6501 0.5653 0.5264 0.4911 0.4755 0.4539 0.4318 
0.7 0.7317 0.7075 0.6189 0.5687 0.5301 0.4897 0.4648 0.4613 0.4055 0.4026 
0.6 0.6995 0.6318 0.5663 0.5195 0.4533 0.4326 0.4019 0.3749 0.3567 0.3367 
0.5 0.6402 0.5842 0.5038 0.4658 0.4286 0.3875 0.3421 0.3191 0.2929 0.2841 
0.4 0.5954 0.5522 0.5028 0.4123 0.3740 0.3565 0.2989 0.2651 0.2437 0.2298 
0.38 0.53870 0.4913 0.4437 0.3856 0.3260 0.2916 0.2620 0.2316 0.2010 0.1837 
0.2 0.5022 0.4578 0.4273 0.3598 0.3321 0.2795 0.2300 0.1747 0.1520 0.1252 
0.1 0.5170 0.4616 0.4006 0.3424 0.3132 0.2611 0.2164 0.1553 0.1106 0.0782 


(a) 
“ 


1.0 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 
1.0 1.0000 0.8684 0.8660 0.6771 0.5656 0.5652 0.5111 0.6327 0.4786 0.5118 
0.9 0.9335 0.8894 0.9131 0.7232 0.5467 0.5089 0.5093 0.5697 0.4969 0.4564 
0.8 0.7918 0.7649 0.6294 0.6942 0.5991 0.4490 0.4575 0.4608 0.3869 0.4027 
0.7 0.7996 0.7652 0.6721 0.5689 0.4886 0.4718 0.4056 0.4858 0.3228 0.3420 
0.6 0.8140 0.7501 0.6625 0.5617 0.5076 0.4206 0.3724 0.3551 0.3605 0.3261 
0.5 0.7978 0.7256 0.5948 0.5663 0.5091 0.3878 0.3282 0.3337 0.2370 0.2504 
0.4 0.8344 0.7665 0.6332 0.5491 0.4925 0.4023 0.3317 0.2549 0.1905 0.2046 
0.38 0.8675 0.8142 0.6973 0.5885 0.4942 0.4097 0.3201 0.2617 0.1907 0.1702 
0.2 0.8592 0.8078 0.6670 0.6029 0.5355 0.4258 0.3402 0.2465 0.1663 0.1066 
0.1 0.8505 0.7491 0.6472 0.5872 0.5286 0.4294 0.3623 0.2548 0.1670 0.0835 


Table 3.27 The numerical values of face area. 


§ 3.13 Voronoi tessellation in higher dimensions 


In § 3.13 it is mentioned that the volume of a tetrahedron is |x;,y;, z;,1|/6, for which the 
absolute value must be taken before adding two or more together. This is what to be found in the 
existing literature. But I think that this definition is flawed because, for one thing, it does not work 
for the case of two dimensions. Imagine what happens if we let a unit cube to always have a unit 
volume. It works as shown in the following. What I discovered here could well be new. 


Assumption 3.4. A unit cube has a unit volume in all dimensions d, d > 2. 


But, apart from the fact that it does not work for the case of two dimensions, the volume 
equation above still appeals to our commonsense, because it looks symmetrical, so we assume further. 


Assumption 3.4. The equation for volume of a d-simplex is V = abs(|x;;,1|)/k, where «;;, 
1 < j < d, is the j:th coordinate of the i‘ point, 1 < i <d+1, and 1 is a unit vector of (d + 1) 
dimensions. 


Under Assumption’s 3.4 and 3.4, and by trial and error, a conclusion is found that the volume 
of a general d-simplex must be such that the & Assumption 3.4 is d!. Therefore we arrive at Definition 
3.1. I actually tried d(d —1) and 24 — 2 before looking at d! and realised that this is the solution for 
k. 


Definition 3.4. The volume of a d-simplex is abs(|x;;,1|)/d!, where xj; is the j*° coordinate of 
the i” vertex of the simplex. 


We can only test the appropriateness of our definition, i.e. Definition 3.4, against the cases of 
two and three dimensions, since four dimensions and above are unfamiliar grounds. But for cubes 
or hypercubes of two dimensions and above, the fit looks encouraging, as can be seen in Table 3.28. 


d 
8 2S Ce OF 8 9 10 
fA Se My ae OM (1) (1) (1) 
£2 4 8 16 32 64 128 (28=256) (29= 512) (210 = 1,024) 
3 9 27 81 243 729 2,187 (38=6,561) (3°? = 19,683) (31° = 59,049) 


Table 3.28 Volume of d-hypercubes whose dimension is £. Values in brackets are implied not 
calculated. 
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The above is, to put it in other words, our attempt at defining the equation for the volume of 
d-simplices and the test of this equation against squares, cubes and d-hypercubes. The equation of 
Definition 3.4 is also rather appealing for the reason that it allows us to calculate the volume of the 
cube as V = abs(|x;;, 1|)/d! = €¢ = [](Az;), where @ is the length of the side of the hypercube and 
the product is over all coordinates xj. 

Table 3.28 justifies our intuitive equation regarding the volume of hypercubes which says that 
such volume is V = £¢. 

The program used in doing the simulations is shown in § A.29. There is also the code for 
testing the formula, which automates the generation of vertices of hypercubes before calculating 


their volume. 
10° : : : : : 10° 


ob 


number of inner cells 
number of inner vertices 


Ll Ll Ll Ll Ll Ll Ll L Ll Ll 
25 3 3.5 4 45 5 2 2.5 3 3.5 4 45 5 


dimension dimension 
(a) (b) 


Figure 3.51 Number of generators required for ten inner cells. The dotted line in (a) has the 
equation y = 3.7 exp(0.98x) while for (b) y = 0.25 exp(2.32). 


The general procedure is to generate the Voronoi cells, exclude boundary vertices, and then 
boundary cells (compare Algorithm 3.3, page 104). This means that producing ten usable cells, 
that is to say, cells in the inner reach, requires for each dimension approximately the number of 
generators plotted in Figure 3.51 (a). Figure 3.51 (b) is the corresponding number of vertices. 

For a Voronoi cell in d dimensions we may draw a box around them such that the walls of the 
box are all parallel to the coordinate planes. For the plane, a Voronoi graph consumes some 60 per 
cent such binding rectangle. For higher dimensions, this ratio goes down from 0.6 to 0.33, 0.15 and 
0.6 respectively for 3, 4 and 5 dimensions. The trend for the standard deviations is also similar to 
this. Both the mean and the standard deviation decrease with dimension, but their curves seemed 
difficult to define if we try to plot it on a semilog graph. The curve turns out to be parabola, as 


shown in Figure 3.52. 
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(a) (b) 
Figure 3.52 Mean and standard deviation of the ratio Voronoi per defining cube. Both graphs 
fit well with the parabolic equation y = ax? + bx + ¢; (a) has a = 0.046, b = —0.498 and c = 
1.418, while (b) has a = 0.00875, b = —.054 and c= 0.187. 


Additional results on a 6-dimensional VT generated from 1,200 generators are p(V,/V,) = 
0.0252 and o(V,/V;,) = 8.5 x 10~*. The results are shown in Table 3.29. 
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Ng Cin Gd BW(Vp/Vn) o(Vp/Vn) Vin 

25 10 2 0.6041 0.0937 33 

30 16 2 0.6034 0.0769 44 

70 13 3 «0.3331 0.0583 241 

200 11 4 0.1488 0.0238 2,374 
500 9 5 0.0638 0.0098 

500 11 5 = 0.0625 0.0130 26,893 
1,200 4 6 0.0252 8.5x10-4 306,710 


Table 3.29 Voronoi volume in higher dimensions; Vp is the Volume of the polytope, V, that of 
the hypercube and ng number of the generators. 
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The accompanying curve plotted in Figure 3.53 is 
a parabola y = 0.03542? — 0.4282 + 1.32. On the 
other hand, neither the mean value nor the stan- 
dard deviation of the volume ratio can possibly fol- 
low a parabolic curve, for the obvious reason that 
both decrease monotonically towards zero. 
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——+ Figure 3.53 The ratio V,/Vp, dimension from 1 
5.5 6 to 6. 
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§ 4. Percolation 

The 1970’s and the 1980’s saw a proliferation of variations on the theme of percolation. Every 
year there seemed to be a new percolation problem or two. For a Bethe lattice Chalupa et al (1979) 
reported a bootstrap percolation where those randomly occupied sites with less than m occupied 
neighbours are recursively emptied one by one until a stable configuration is reached. The problem 
they are interested in is that where the impurity concentration, dilution and crystal-field interaction 
compete in magnetic materials compete against the exchange interaction, resulting in the magnetic 
moments and consequently the magnetic order being destroyed. 

Wilkinson and Willemsen (1983) introduced invasion percolation. Working with Schlumberger 
they were interested in the real problem in the oil industry where water displaces oil via capillary 
action. Their approach was that of a constant flow rate, not one of a constant pressure as usual. 
Here water displaces oil in the smallest available pores. But when water completely surrounds any 
region of oil, no further advance into that region will be possible, oil being incompressible. Such 
regions are called trappings and cause a problem generally known by the name of residual oil, an 
economic bane for oil industry. 

In the above example the hydrophobic versus hydrophilic propertyplays an important role in 
the replacement of oil in pores with water. And water is prevented from penetrating trapped oil 
regions by much the same principle as that which prevents the water in the contents of a sandwich 
from crossing the spread layer of butter to wet the hydrophilic bread. For this latter case the pores in 
question are those within the bread texture, and the soaking of the bread is best prevented provided 
that the trapped regions of butter or margarine percolates in two dimensions to form a single layer 
or film which entirely covers the sectioning surface of the bread. Moreover, there is a similarity 
also in the internal structure of both the rock and the bread. The density of pores in bread is not 
homogeneous as a result of the tension on the surface of the dough caused by internal air pressure 
originated from the yeasts inside, as well as because of the heat applied to it when inside the oven, 
which dehydrolises the surface and makes it dry and hardened. The same inhomogeneity can be 
found inside the rocks which form the oil reservoirs where the regions of oil are surrounded by rocks 
which are less porous and have less permeability. 

For Adler and Aharony (1988) a random walker, aka ant, treads on clusters. The ant enlarges 
a cluster by stepping on to an empty site next to it which meets certain conditions. They called 
this problem diffusion percolation. An example of a condition met is where the empty site has two 
or more occupied neighbours on a square lattice. 

Most percolation in studies happens by randomly toggling the phases of sites or bonds in regular 
lattices. Kerstein (1983) considered randomly located spheres, take the complementary region of 
their union, and then perform percolation on the former. He showed that such percolation problem 
is equivalent to a percolation on a Voronoi lattice whose sites are the sphere centres. 

In the same way as an infinite loop in computer science means that one can not come to the 
termination of a program going along the time dimension, an infinite cluster in percolation means 
that one can not come to the end of a cluster shifting along either one of the dimensions. The 
former case is along the one dimension of time and is only possible because of the time flows in 
one direction. Therefore half the line spanned is considered as being infinite. In one nondirectional 
dimension, except for the trivial case where one can consider the entire space as being one single 
cluster, no infinite cluster is possible. In simulation, when a cluster spans the whole of the space 
considered along any dimension we say that it is infinite since as one moves a cross section along 
that dimension, it always contains a section of the cluster. 

The bond percolation program of two-dimensional Voronoi network by Tiyapan (1995, KNT(iii), 
p. 78) takes O(n?) time provided we assume that the contribution made by the number of clusters 
together with the that by each cluster comparison amount to a linear time complexity term. In order 
to justify this assumption, consider first the number of clusters. The maximum number of clusters 
possible depends on the size of the system, in other words it must vary as nf, where 0 < k, < 1. 
This maximum value should be approximately 0.5 because of symmetry between occupancy and 
void clusters. Next consider the time involved in comparing two clusters. On Matlab this is a sparse 
vector comparison which is likely to involves some linked lists, and similarly should take time as a 
function of nk, 0 < ko < 1. Because on average the size of clusters is always small before Pc, ke 
will be less than 0.5. Therefore n¥ - nk < n and it is safe to assume O(n) time from both of them 
combined. g.e.d.. A C translation (Tiyapan 1995, ibid., p. 80) of this program, though not as bad 
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as it may seem because constant coefficients are small, gives O(n®) time in comparison. 

In the field of geology Miller et al (2000) studies the analogy between the dilatant slip in 
earthquakes and the hydrofraction occurred in melting and dehydration, the percolation of the 
latter in the permeability network internal to the fault being the cause of the former. According 
to them, the existence of toggle switches in the permeability rules out the assumption that the 
permeability throughout the whole system is homogeneous. Instead, the system reorganises itself 
into sytems of different scales of interaction according to the degree and nature of its inhomogeneity. 
At the critical state the scale of interaction is equal to the scale of the model. 

The percolation probabilities given by Stauffer and Aharony (1994) are, the first number being 
for sites and the second for bonds, for the honey comb lattice 0.6962, 0.65271; square 0.592746, 
0.50000; triangular 0.50000, 0.34729; diamond 0.43, 0.388; simple cubic 0.3116, 0.2488; body-centred 
cubic 0.246, 0.1803; face-centred cubic 0.198, 0.119; 4-d hypercubic 0.197, 0.1601; 5-d hypercubic 
0.141, 0.1182; 6-d hypercubic 0.107, 0.0942; 7-d hypercubic 0.089, 0.0787. 

De Gennes et al (1959) investigated disordered binary solid solution AB where active atoms A 
are randomly replace the nodes of the periodic matrix B. There exists a critical concentration A in 
B below which all clusters are finite, and above which both finite and non-finite, i.e. infinite, clusters 
exist. Such solid solution in networks can represent the spin waves in alloys with one ferromagnetic 
component or the impurity bands in semiconductors. They cited seminal work on percolation by 
Broadbent et al (1957), but no mention was made about the Ising model. 

In a way, percolation is similar to diffusion. In diffusion the particles considered move about 
randomly, whereas in percolation they can only crop up randomly at predetermined locations on a 
network which is fixed. We could imagine, for instance, cars running along the roads within a traffic 
network as diffusing through them. Then the percolation could occur on a larger scale, that is the 
scale of aroad. The cars move along, that means they diffuse; but the roads remain fixed, and so their 
phases could percolate. In other words, in diffusion the particles move while in percolation, whether 
there are moving particles or not, it is the phases that percolate. Since historically percolation began 
as the study of diffusion of particles in a network of tubes in which the phases are naturally defined 
as the tubes being blocked or unblocked, these definitions have become most frequently used in other 
areas of application, for example in filtering membranes and traffic networks. 

But this is not necessarily the case. Instead of dealing with a fixed network, one may consider a 
model of percolation in a continuum, for example by randomly patching an area until all the patches 
connect with one another somehow and percolate. The patch could be of any shape, as well as 
polygonal and circular. We can consider the percolation in a certain area as having occurred when 
there appears a cluster of patches which traverses any two opposite sides. One application of this is in 
the study of occurrences of epidemics. Hoyle and Wickramasinghe (1979), having given a convincing 
argument in favour of viruses and various forms of diseases being carried to Earth from space 
by comets, talk about patchiness of pathogenic clouds. According to them, simultaneous attacks 
across vast region rules out person-to-person theory. Moreover, influenza epidemics are generally 
characterised by sudden onsets and equally sudden ends. These epidemics and plagues may be 
thought of as the percolation of these patches in a sufficiently large, predefined area. The bacteria 
and viruses coming from space adding to gene give the possibility of jump patterns in evolution (cf 
Smith and Szathmary, 1995). The cells deliberately refuse to block viruses because they could prove 
to be useful in the long run, generally not by individuals but by the species. Historical examples are 
the disease described by Thucydides between 431 and 404 BC, five epidemics of ‘English Sweats’ 
between 1485 and 1552, and more than ten influenza pandemic from 1700 to 1900. Example of 
diseases which are caused by bacteria and viruses are bubonic plague, chicken pox (varicella), cholera, 
common cold, Legionnaires’ disease, leprosy, measles, mumps, poliomyelitis, small pox, tuberculosis, 
and trachoma. Examples of major evolutionary transitions are those going from RNA to DNA, from 
prokaryotes to eukaryotes, and from asexual cloning to sexual propagation. Another example is the 
transition from primate to human both of whom differ from each other neurologically in the ability 
to use language and the power to conceptualise. One description of the Great Plague in London 
(Dickens, 1851). There was a rumour that a few people died in the winter of 1664. In May 1665 the 
disease burst out in St. Giles’s which raged through July and September in every part of England; 
approximately 10,000 people died in London alone. But then the equinox winds virtually blew the 
disease away, and the Plague quickly disappeared. The existence of interstellar organic matters is 
supported by strong evidences with more and more complex substances constantly found (cf Hoyle 
and Wickramasinghe, 1978). 

Gas turns into liquid by the growth of larger and larger clusters. But unlike percolation in 


Ph.D. Thesis, UMIST. K N Tiyapan. Chapter 4: Percolation 


networks, clusters in condensation process are not well defined. Crystalisation is also characterised 
by the growth of one phase within another. But here the orientations are an inherent part of the 
clusters themselves, and there is a long range coordination among clusters. One important thing in 
percolation is for the system size to be infinite. Random discs overlapping one another in a continuum 
helps correct counts of bacteria cultures (cf Essam, 1980). Let p, be the critical probability and 
P(p) the percolation probability. Then as p approaches p, from above, P(p) & (p — p.)° where 
0.4 < 8 < 0.5 is the critical exponent. If p approaches p, from below, S(p) & (pe — p)~7 and 
&(p) © (pe —p)~” where S(p) is the conditional average (s|F'), s the number of particles in a cluster, 
and F the event that the origin O is occupied and belongs to a finite cluster. In other words, 
S(p) is the mean size of cluster at the origin given that F’ occurs. This last event occurs with 
probability pr = p(1 — P(p)) and pr = p for p < pe. Furthermore, S(p) = Dr >, C(r, p) and 
£?(p) = [pr S(p)]~! 4, r?C(r, p) where the pair-connectedness function is C(r,p) = (n(r)|F)pr, and 
n(r) is the indicator defined to be one if r is connected to O and zero otherwise. When F occurs, 
s = )0,n(r). There are two different definitions of the critical probability, p, = sup{p|P(p) = 0} 
and 7. = sup{p|P(p) = 0 and S(p) < oo}, sometimes denoted by py and pr for Hammersley and 
Temperley respectively. By definition, 7, <p... When p approaches p, from above, 7’ and v’ are 
similarly defined respectively by S(p) * (p — pc)” and €(p) © (p—p-)” It is generally assumed 
as y' = y and v' = v. No proofs for these assumptions exist, even though they are consistent 
with the series expansions. Estimates of y and v are 1.6 < y < 1.7 and 0.8 < vy < 0.9. Ina 
dilute ferromagnet a cluster containing s spins each of which has a unit magnetic moment has a 
probability p = exp($sh)/ [exp(—4sh) + exp($sh)] of being parallel to the magnetic field H > 0. 
Here h = 2H/kgT where kgT is the thermal energy. For an infinite cluster, p = 1. The zero-field 
magnetic moment is po(p) «x P(p). The zero-field magnetic moment is po(p) «x P(p) ~ (p—p-)8. The 
field dependent magnetic moment at p- is ((1— exp(—sh))/(1 + exp(—sh)))r and $(1—G(pe,h)) < 
pc(h) < 1— G(pe, h) where G(p, h) = (exp(—sh))”. Therefore the critical probability 6 is defined 
by pic(h) ~ 1— G(pe,h)h'/°. The correlation function between o on site i and j is defined by 
Ti; = (o10;)7 — (0%)r(o;)r where (-)7 is an average over spin states. Then Tj; = 1 when i and j 
belong to the same cluster and zero otherwise. The fluctuation formula is kpT (x) = >); (Ti) where 
(ij) = C(r; — ri, p), the mean susceptibility is (x) = prS(p)/kaT ~ (pe — p)~7 and the mean free 
energy is (F) = (kgT In2)K (p) where K(p) is the mean number of clusters and K (p) ~ (p. — p)?~°, 
where the index assigned to the third derivative divergence is 1+ a, a % —0.5. 


The growth mechanism of clusters in percolation is all so found in biology. Williams and 
Bjerknes (1972) simulate a tumour in the basal layer of an epithelium. The basal cells become less 
sensitive to Charlone which controls cellular division, and thus they divide « times faster than the 
normal cells where & > 1 is the carcinogenic advantage. Abnormal cells interior to the basal layer 
divide, push, and then replace the neighbouring cells leaving the overall configuration unchanged 
except at the border where the abnormal cells exert a thrust of « — 1 on their normal neighbours, 
that is av = (« —1)n where n and WN are respectively the numbers of peripheral and total abnormal 
cells. They found that dimensionality of fractal is involved and the dimension 1.1, instead of 1, must 
be assigned to the periphery, which means that n is proportional to N°> not N°5. They found 
that abnormal cells push out faster in this order: the triangular, square, and hexagonal lattice. We 
may explain this, by looking at the coordination numbers of these three lattices, that the higher 
coordination number the lattice sites have, the slower the cluster expands. Added coordination 
means a higher degree-of-freedom the newly divided cells have to move about while still remaining 
local. 


Percolation is a field full of far more open questions than discovered answers. This is one of 
the reasons it is more suitable to research than university examinations (cf Stauffer and Aharony, 
1998). The fascination not so well hidden within the field is reflected in the number of research 
monographs written on this field that is close to percolating, which makes it impossible to give a full 
list of references. Percolation deals with the clusters formed by randomly occupying the each site 
of a very large lattice with probability p independent of its neighbours. Most of the studies to date 
concentrate on the critical phenomena, which exist in a very narrow range in each problem, and the 
scaling theory which tries to describe them. The contiguity criteria are easier to define in applications 
in physics and physical sciences, because here problems can be described by, or easily simplified into 
some definite geometry, whereas in other applications this may not be the case. For example, in 
economic modelling the nature of the geometry or even the number of dimensions of such structure 
is still not well understood (cf Tiyapan, 1997, KNTs(vii) and KNTs(viii)), let alone the problem of 
contiguousness. Traffic network (cf Tiyapan, 1997, KNTs(i) and KNTs(ix)) is another example. Here 


113 


114 Ph.D. Thesis, UMIST. K N Tiyapan. Chapter 4: Percolation 


the flow is along channels within tubes, i.e. roads, both of which could be directed, which make 
the network a pool of tangled threads not only difficult to imagine but also to define such things as 
congestion and contiguousness. Even in geographical problems like that of forest fire, the spread of 
fire may be defined in at least three ways, namely by next-nearest neighbour, neighbour or double 
neighbour respectively in cases of a common corner, a common side or two neighbours. Reservoirs 
occur when the petroleum, formed in sedimentary rocks, migrates into permeable sedimentary rocks 
like sandstone. Around 70% of North Sea oil fields also date from the Jurassic period. These include 
Beatrice, Brent, Cormorant, Murchison and Tartan. Oil companies want to exploit reservoirs where 
p > pe. Moreover, they want to tap into the largest cluster in these reservoirs. For this purpose 
bore hole samplings are carried out which collected within L x L frames and then the number of 
points, m(L), belonging to the same cluster counted. It turns out that m(L) is proportional to 
L'-®. In general m « L4 where d is either integral, in which case it is the Euclidean dimension, or 
nonintegral, in which case it is the fractal dimension introduced by Benoit Mandelbrot. There exists 
a correlation length € such that m(L) « L*® for L < €, and m(L) « L? for L > €. This correlation 
length is a measure of the largest hole of the largest cluster and decreases as p is increased above pe. 


One way of finding the percolation threshold is by using the ant in a labyrinth algorithm shown 
here as Algorithm 4.1. 


Algorithm 4.1 Ant in a labyrinth. At p., k = 1/3, whereas k = 0 for a constant distance, and 
k =1/2 for normal diffusion. 


for various p do 
repeat a large amount of times 
occupy the sites with probability p; 
sum < 0; 
for a large number of simulations n do 
for various ¢ up to a large number do 
sumr + 0; 
fori=1tot 
ant randomly occupies an occupied site; 
identify neighbours; 
move into an occupied neighbour randomly chosen; 
r < distance traveled; 
sumr <— sumr +73 
endfor 
sumrsq + sumr?; 
endfor 
endfor 
sum <— sum + sumrsq; 
R¢sum\/?; 
plot R against ¢ in double logarithmic scales; 
endrepeat 
k + the slope of the straight line just plotted; 
endfor 
De < the p such that k = 1/35) o 


Also, the number of steps the ant takes, ¢, for the linear size of a region it visited, R, is 
fractal, that is t x R!/*. They give several percolation probabilities for selected lattices, which are 
shown in Table 4.1, as well as the exact values of p,’s, which are for the square bond percolation 
1/2, triangular site 1/2, triangular bond 2sin(z/18), and honeycomb bond 1 — 2sin(7/18). For the 
honeycomb site percolation, p. < 1/2. Site percolation on hypercubic lattices of high dimensions 
have p, = 1/(2d— 1). 
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lattice site bond 
honeycomb 0.6962 0.65271 
square 0.592746 0.50000 
triangular 0.500000 0.34729 
diamond 0.43 0.388 
simple cubic 0.3116 0.2488 
body-centred cubic 0.246 0.1803 
face-centred cubic 0.198 0.119 
4-d hypercubic 0.197 0.1601 
5-d hypercubic 0.141 0.1182 
6-d hypercubic 0.107 0.0942 
7-d hypercubic 0.089 0.0787 


Table 4.1 Percolation thresholds from Stauffer and Aharony (1998) 


There is no percolation in one dimension because in such case the percolating cluster must 
necessarily contain the whole space. But there are some applications where the critical blockage is 
important, for example the blockage of drainage grilles by pea shingle is one-dimensional. Here the 
critical amount of blockage depends on the critical rate of flow which in turn depends on the amount 
of water and the rate of accumulation of water to be drained. The maintenance of gullies and grilles 
is done by cleaning, flushing and grit bucket emptying (cf Harrison and Trotman, 2002). 

For the Bethe lattice p, = p. = 1/(z—1). In general, p. < p,. For a network of d dimensions, 
zpe & d/(d—1). In his review, Sahimi (1994) gave some of the current know p,’s which are listed in 
Table 4.2 here. 


d Z De 2ZPe iP 

honeycomb 2 3 1 -sin(7/18) 1.96 0.6962 
square 2 4 = 1/2 2 0.5927 
kagome 2 4 0.522 2.088 0.652 
triangular 2 6 2sin(7/18) 2.084 1/2 
diamond 3 «4 0.3886 1.55 0.4299 
simple cubic 3. («C6 0.2488 1.49 0.3116 
body-centred cubic 3 8 0.1795 1.44 0.2464 
face-centred cubic 3 12 0.198 1.43 0.119 


Table 4.2 Percolation probabilities, cf Sahimi (1994). 


The accessible fraction, fy, is the fraction of occupied bonds that belong to the infinite cluster. 
The backbone fraction, fg, is the fraction of those accessible bonds which belong to a transport path, 
i.e. a path with all dead ends excluded. The correlation length, A, is the length scale over which 
the random network is macroscopically homogeneous. In a Monte Carlo simulation it is necessary 
that the size L of the network is sufficiently larger than this correlation length in order to obtain 
a pe which is independent of L. The expected cluster size m is E(m) = 2, m?mm/Qom, Mm); 
where n(p) is the expected number of clusters of size m per lattice site and mn,, the probability 
that a site belongs to an m-cluster. The fraction of the network which can accommodate flow has 
various names associated to each application, for example the effective electrical conductivity ge, 
the effective diffusivity D, and the hydrodynamic permeability k. The effective elastic moduli, G, 
are the elastic moduli of the network a fraction p of bonds of which are elastic elements, while the 
rest are rigid or stiff. The fraction of isolated occupied sites is f;(p) = p — fa(p). The universal 
scaling laws state that Pe(p) ae (p — pe)", fa(p) oe (p— pe)°?, fa(p) ae (p—pe)**, A(p) id |p — Pe|~”?, 
k(p) ~ |p — pel~”, ge(p) ~ (p—pe)# and G(p) ~ (p—p-)4, where bg, bp, Vp and 7p are topological 
exponents, completely universal, depend only on the dimensionality but not the microscopic details 
of the system. The effective diffusivity D.(p) ~ (p — pc)#—"? because ge ~ NeD, and ne ~ fa(p). 
Near pe, k(p) ~ (p— p-)®. For network percolations e = y, but for the continuum percolation this 
may not be the case. 


As p > pz, perfectly conductive clusters become larger and g, increases. Then g.(p) ~ 
(pe — p)~* near pe and diverges at p-. In two dimensions, 4 = m. If p edges are totally rigid while 
the rest are elastic, then G diverges as p + p> such that G ~ (p, — p)~$ For large m near pz, 
Nm ~ m-" f [(p— pe)m??], where 7, and op are universal and f(0) 4 0. Some relations among the 
geometrical exponents are tT) = 2+ bpop and vpd = bp +1/op = 2bp + yp. Sahimi (ibid.) listed some 
values of the current critical exponents and fractal dimensions which we list again here as Table 4.3. 
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d=2 d=3 Bethe lattices 


bp 5/26 0.41 1 
bg 0.48 1.05 2 
Up 4/3 0.88 1/2 
‘p 43/18 1.82 1 
op 36/91 0.45 1/2 
Tp 187/91 2.18 5/2 
D. 91/48 2.52 4 
Dg ‘1.64 1.8 2 
Dig: WAZ 1.34 2 
bu 1.3 2.0 3 
m 1.3 0.73 0 


Table 4.3 Critical exponents and fractal dimensions, cf Sahimi (1994). 


The system is macroscopically homogeneous when L >> A but not when L < X where the cluster 
spanning the sample is self-similar and fractal at all length scales up to A. Its mass is M ~ \”-, 
where D, = d— bp/vp is the fractal dimension of the cluster. For L > », D. = d. Similarly the 
backbone is also a fractal object when L < 4, and it has the fractal dimension Dg = d— bg/vp. 
When L < ) one should replace \ by L, at diverges at p= p- and the red bonds have M, ~ L?*, 
D, =1/vp, when L < X. The minimal or chemical path between two points of a percolation cluster 
is the shortest path between the two points, and Lmin ~ L?™™ for L < x. 

Finite-size scaling is the effect of finite size on the critical properties. As p—> p., AX > L, and 
the variations of a property P_ becomes Py ~ L~*f(u), where u = L1/¥»(p —p,) ~ (L/A)/” and 
f(0) #0. When p > p. and L + 00, Px ~ (p—pe)® and « = 6/vp. There is also a shift in the 
percolation threshold, and p.—p,(L) ~ L~!/” where the percolation probability is p, for the infinite 
system while its effective value for a finite system is p-(£). Correction should be made whenever 
the system simulated is small, that is Pp, ~ L~* [a; + aggi(L) + aggo(L)], where g; and go are the 
correction terms to the scaling. For transport properties, for example conductivity, diffusivity and 
elastic moduli, g, = (In Z)~+ and go = L~! are often good enough. 

Sahimi (1994) mentions a random distribution of inclusions, such as circles, spheres, or ellipses, 
in an otherwise uniform system or continua; a percolation on random polyhedra such as the VT; and 
also a random distribution of random conducting sticks of with given aspect ratio, or plates a given 
extent, which is important in modelling fracture networks in rocks. In the study of percolation in 
continua, one considers a continuous, random function h(r), defined for all points r’s in the entire 
space, such that (h(r)) = 0. All regions of space where h(r) < R are in a phase which is different 
from the rest, the volume of which goes from zero to infinite as R changes from negative- to positive 
infinite. As R becomes larger, so do the islands, and the phase grows and percolates at the critical 
R,.. For continua percolation, the critical occupied volume fraction is defined as ¢, = p.f, where f 
is the filling factor of a lattice when all its sites are occupied by an impermeable sphere in such a 
manner that those which are nearest neighbours touch each other. This volume fraction seems to 
be an invariant, approximately 0.45 for two- and 0.15—0.17 for three dimensions. If the spheres are 
permeable, then ¢, = 1— exp(—B,./8) where B, = zp, is the number of bonds per sites at p., and 
in continua percolation B, > p.z as z — oo. If ¢, is an invariant of continuous systems, then R, is 
given by ¢, = Wie h(r)dV, where V is the volume of the system. 


§ 4.1 S-curves and the percolative phenomena 


Before any abrupt change occurs there must be a graduation process leading up to it. In 
processes where the time constant is long, the characteristic s-shape is obvious, take for example 
a learning curve. But even where time constants are short, it is doubtful whether such thing as a 
strictly abrupt change exists. At the very point of transition, undoubtedly there may be a singularity, 
for instance at the Big Bang. But even the Big Bang can not exist alone by itself. To preserve the 
symmetry of nature, there must be another process at the other end leading up to it. It is only 
because such process must necessarily be on the other side from us and we can not see it from 
here. That side may belong to the antimatter or the anti-universe, but I believe whenever there is 
a singularity there must be a symmetry. 

An earlier work by Tiyapan (2003, KNT8(iii)) can be used as an example. The per cent extraction 
curves which I drew then (Tiyapan, 1991, KNT1(i)) all show the change to be abrupt, starting off 
from zero time immediately with a positive gradient. This can not happen in the real world, so 
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there must have been a foot of the characteristic s-curve at the beginning. It must have been that 
the sampling time used is too long, and anyhow the nature of the reaction may make it impossible 


to observe the development in detail with accuracy. 


0.97 
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monoatomic solids according to Debye. 


The foot of a positive s-curve has a positive, nonzero second derivative. This corresponds 
with the positive cooperativity of the product curve in studies of enzyme assay and kinetics. In 
enzyme assay, the product versus time graph shows a positive cooperativity characteristic when 
the Hill constant h > 1 in the equation y, = [s]"/(K + [s]"), where y, is the fractional saturation 
of the enzyme with substrate while s the concentration of the substrate (Eisenthal and Danson, 
2002). Cooperativity reflects the equilibrium binding of substrate or other ligand. The binding of 
a substrate molecule to an enzyme either facilitates, when the cooperativity is positive, or hinders, 
when the same is negative, the binding among molecules of the same substrate. 


§ 4.2 Voronoi percolation in two dimensions 
The percolation of Voronoi tessellation in two dimensions can be achieved by the following algorithm 


Algorithm 4.2 Network percolation in 2-d. 


generate random points; 
generate Voronoi tessellation and Delaunay triangulation; 


find vertices within the square box bounded by lower and upper bounds; 
find neighbours of all cells, bonds, vertices, and edges; 
for unit = cell, bond, vertice and edge do 
find number of units; 
permute list of units; 
clear cluster list 1, cluster list 2, set of resultant clusters; 
cluster percolated < false; 
for i= 1 to number of units do 
existing cluster joined = false; 
for 7 = 1 to number of clusters in cluster list 1 do 
if the j* cluster contains the i” unit in permuted list do 
merge the i” unit into the j*" cluster ; 
existing cluster joined < true; 
end 


if existing cluster joined is true 
move the j* cluster of cluster list 1 to cluster list 2 ; 


for k = 1 to number of clusters in cluster list 1 do 
if the k*” cluster in cluster list 1 touches the cluster in cluster list 2 do 
merge the former into the latter; 
else 
append the former to cluster list 2; 


end 


end 
test percolation of the cluster just updated; 


move cluster list 2 to cluster list 1; 
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clear cluster list 2; 
break; 
end 
end 
if existing cluster joined is false 
create a new Cluster of size one and append it to cluster list 1; 
end 
append cluster list 1 to set of resultant clusters ; 
end 
end 


The Pc probabilities are found to be (2) “4/2 Pc = 0.511040.0856, — 77,"p» = 0.3095+40.0523, 
Vn 2dp, = 0.7231 40.0616 and ¥",24p, = 0.6801 + 0.0468. 


And the coordination numbers are —.2436 7, "4@c = 5.2320, — .2979 45, *4a5 = 9.5022 , 
0259 ¥" *4n, = 2.8617 and 0382 ¥" 4x, = 3.8064. 


The percolation theory in two dimensions has benefited much from discoveries regarding the 
Ising model in Physics. The Ising model is to date probably the most successful model in percolation 
theory. This is not only because it can describe in details the phenomena of ferromagnetism and 
antiferromagnetism, but also because it does so by replacing a noncrystalline solid in the atomic 
scale with a perfect lattice. 

When a bar of iron is placed in a external magnetic field at a constant temperature, the field 
will induce some magnetisation in the bar. If the external field is slowly turned off, there are two 
scenarios possible. If T < T, the bar retains some of its internal magnetisation, but if T’ > T, the 
magnetisation completely disappears. The transition temperature T, is known as the Curie point. 

The Ising model represents the iron bar by lattices, for instance a square lattice. Then it 
defines two power series, namely the low- and the high temperature series. The low temperature 
series comes from a bivariate generating function which is generated from the number of colourings 
for which there are p black sites and q black-white edges. In other words if we let A(p, q) be the said 
number of colourings, then our bivariate generating function becomes a(z,y) = pa é A(p, q)uPy4 
and the power series is obtained from a(z, y) limy—o.(1/N) In(a(a, z)). 

On the other hand, the high temperature series of the Ising model changes its form with the 
dimension. In one dimension it degenerates, while in two dimensions it is isomorphic to the low 
temperature case and comes from a bivariate generating function that is generated from the number 
of even polygonal drawings whose area is p and which has q edges. In three dimensions it comes from 
a univariate function that in turn is generated by the number of even polygonal drawings which have 
q edges. An even polygonal drawing on a lattice is a union of its subgraphs that uses each edge of the 
latter at most once and each site an even number of times. It is also known as an Eulerian subgraph 
and is indeed a union of simple, closed and edge-disjoint polygons which need not be connected. 

In three dimensions if we let B(q) be the number of even polygonal drawings mentioned, then 
our univariate generating function becomes b(z) = )/, B(q)z? and the power series can then be 
obtained from f(z) = limy_,.(1/N) In(b(z)). 

Two main problems remain unsolved regarding the Ising model, namely that of finding closed- 
form expressions, of a(x,y) for two dimensions and of 6(z) for three. 


§ 4.3 Voronoi percolation in three dimensions 


In three dimensions the procedure of finding p, is similar to that used in the 2-d case. Algorithm 
4.3 written was used in doing the simulation. Developed for running on Matlab, it does every thing 
in matrix form. As a consequence, there are various types of data all of which are matrices. These 
data structures can be summarised into the following. 

A list is a one-dimensional matrix whose dimension is the number of its members. An indez is 
a matrix whose members are either one or zero. These index matrices in two dimensions are like a 
map or an array. They are normally one or two dimensions, and map the relationship between the 
members of one set and another. A set is an m x n matrix every one of the members a,j; of which 
is a matrix. A pair matriz is a square diagonal index matrix that maps among members within the 
same set. It contains the information about relationships like neighbourhood, group membership, 
connection, etc. A pair matrix or an index can also contain information other than the existence or 
membership flags, for example the edge length in the case of a vertices pair edge matrix. In a bond 
or an edge the mid point represents its position. 
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The program can be divided into three parts, creating and arranging the Voronoi data, finding 
the neighbourhood matrix and the percolation simulation for p,. Putting the vertices of a face in 
order amounts to juggling from one end of each stick, i.e. edge, to another. Edges are traced in one 
direction only until all the edges are successfully linked head to tail. Two lists receive the result 
from the tracing, vertices which match go to one of them while those which do not is put in the 
other and recycled. 


Algorithm 4.3 Managing the Voronoi data in three dimensions. 


xz + create random points; 
(Va, Ca) < create Voronoi tessellation; 
t + create Delaunay tessellation; 
find v, such that for all vg, 0 < ug <1; 
find c, such that all vertices of cg are in Un} 
CH Cy} 
for all do 
mp < find mid bond coordinates; 
b; + find bond length; 
endfor 
for all pairs of cells do 
fa < find shared vertices when either of the two cells is in cy; 
endfor 
for all vz, v — Va3 
for all c,, f <— fa3 
for all f’s do 
if it has three vertices then 
all the three possible pair combinations are neighbours; 
else 
(a,b,c) + find the face equation from three vertices; 
6; + ka where k =1/(a? +b? +c?) and i= 1, 2 and 3; 
max + 0; 
for j = 1to3 do 
if 6 < max then 
jm — 53 
max  6(jm)5 
endif 
endfor 
if j, = 1 then 
pos a7; 
elseif j,, = 2 then 
pt a3 qt 23 
else 
pea oo 
endif 
d+ find delaunay triangulation from p and q; 
for all edges of all triangles of d do 
record the number of times they occur; 
endfor 
neighbours ¢ vertices of edges which occur only once; 
for all f’s do 
order all their vertices; 
end 
endif 
endfor oO 


The terms vertices and edges refer to the Voronoi tessellation only. The vertices and edges of 
the Delaunay tessellation are cells and bonds of the VT. The program perco3d.m can find p, for all 
of these. For each one of them we must find the neighbourhood matrix as well as lists of the upper- 
and the lower boundaries. The order of blockages is completely decided in advance by finding a 
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permuted list of numbered items. The program finds the history of clusters for the whole range of 
p, i.e. from 0 to 1. 

To obtain the matrix of cell neighbours, first crop out between 5 and 10 per cent the boundary 
in all directions. Then the neighbour matrix is found from the DT matrix. 

The bond neighbour matrix is simply the rearrangement of the cell neighbour matrix obtained. 
We have already found the vertices neighbour matrix earlier while searching for vertices shared 
between cells. This is rearranged for the use in the percolation program, and then again for the edge 
neighbour matrix. 

Next is Algorithm 4.4 which carries out the percolation simulation. The cluster information 
is swapped alternately between three variables, A, B and tmp, to facilitate the flow; A(i) is the i” 
cluster in A. All variables are assumed to be cleared at the start of the algorithm. At the end of 
each run the value of p, is computed, and the list p is reversed and reused for a second time if it 
was the first run. 


Algorithm 4.4 Voronoi percolation in three dimensions. 


for each permuted item p(i) in the list do 
joined <+ 0; 
for each cluster A(j) in A do 
if A(j) contains p(i) then 
A(j) © AG) N pla); 
B(1) « AG); 
tmp + A; 
A+< tmp — A(j); 
for all A(k) in A do 
if A(k) A B(1) then 
B(1) + A(k) N BCA); 
else 
B(++np) + A(k); 
endif 
endfor 
percolated ?; 
ACB; 
break; 
endif 
endfor 
if joined then 
A(+ + na) + pli); 


endif 
endfor Oo 
The resulted Pc probabilities from simulation are ("pe = 0.2340 + 0.0448, 1, 84.4 po 
=0.1178+£0.0271,  ¥" 94... p, = 0.2941 40.0831 and "94, = 0.4311 + 0.0324. 


And the coordination numbers are — 3329 43, °4a- = 11.1231, 5996 (3, 34, x) = 23.0548 , 
0245 13 34.) 2v = 3.6926 and 0432 43, °*r. = 5.5002 . 

Four different types of percolation have been performed, namely cell, bond, vertice and edge 
percolation. Vertices are inherently zero dimension, edges and bonds one dimension and cells, which 
in reality have three dimensions, is considered for this purpose as having none. When considering 
vertice- and cell percolation, the position of the vertices and nuclei determines whether they could 
connect somewhere to somewhere else. When the object the percolation of which we consider has a 
nonzero dimension and the dimension of the network is small, that is to say, in the same order as 
that of itself, p. obtained will depend on the network dimension. In other words, there will be an 
influence from the size of the objects on, and distort, the space in which the percolation occurs. 

Additional results on a larger network are, the Pc probabilities from simulation are 7,)'*pe 


= 0.2039 + 0.0410,  ¥",94p, = 0.0963 + 0.0114, ¥7,°4p, = 0.1659 0.0571 and -¥n,84p, = 


0.4172 + 0.0242, and the coordination numbers are v3, °4n_ = 12.4177 , Mee *a5 = 25.7077 
5 ay ay = 3.8126 and =, 4a. = 5.7050 . In this case the network was originally built from 


1,000 Poisson point generators, and for the simulation we have n. = 723, ny = 4,489, ny = 4,855 
and ne = 9,255. 
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§ 4.4 Percolation of 2-dimensional Voronoi sections 


The program in § A.4 does a 2-d section of the 3-d VT. Algorithm 4.5 describes how it works. 
It assumes that the Voronoi tessellation already exists. Here d,, means a denominator while n, a 
numerator, V and C’ means vertices and cells of the 3-d Voronoi tessellation while v and c those of 
the 2-d section. In particular, c € C 


Algorithm 4.5 Plane section of Voronoi in three dimensions 


(Ag, Ay, Az); — (@2 — £1, ye — yi, 22 — 21); for all edges 4; 
for all edges 7 do 
dm + aa + bAy + cAz; 
if d,, nonzero then 
Np — ax, + byi + 6213 
t¢ —ny/dm} 
if0<t<1 then 
(x,y, 2) = (a1 + Az, y + Ay, a1 + Az)a5 
{us} = (2, y,2)5 
endif 
else 
atcate; 
dm + aa + bAy + cAz; 
Np — ax, + byi + €213 
t¢ —n-/dm} 
if0<t<1 then 
(x,y, 2) — (a1 + Az, yy + Ay, a1 + Az)as 
{us} — (2, y,2)5 
endif 
if n, = 0 then 
{vs} & (21,41, 21)3 
{us} - (x2, y2, 22)5 
endif 
endif 
endfor 
for all c, do 
find the Delaunay triangulation; 
count n, of all the triangles, add the numbers into a single list; 
endfor o 


The intersection of a plane az + by + cz +d =0 by the line which is defined by # = x, + Aat, 
y = yr + Ayt and z = 2 + Azt, where Az, Ay and Az are respectively (a2 — 21), (y2 — yi) and 
(zg — 21), is determined by t = —(ax1 + by1 + cz1 + d)/(aiAa + bAy + cAz). If the denominator is 
zero the line is parallel to the plane, and if the nominator is also zero contained therein. 

The results from the percolation simulation on the section V3 are y),95)°pe = 0.5494£0.1223, 
2S py = 0.3515 + 0.0764, 022%, = 0.7557 + 0.0757, Hat”. = 0.6210 + 0.0665. 

The coordination numbers are —_.2212 ate we = 4.6894 , .6021 yi ep = 8.8100 , 
0387 (5 ay = 2.7495, 1118 8, "*e, = 3.6691 . 


When sectioned by the plane (a,b,c,d) = (0.01, 0.5, 0.5, —0.5), our 
VT gives a picture as shown in Figure 4.2. Here the points shown are 
merely the average of the coordinates of all vertices in each cell. From ten 
simulations of these 118 cells, 300 bonds, 288 vertices and 405 edges having 
respectively the coordinate numbers of 5.0847, 9.7933, 2.8125 and 3.7333, 
we obtain p, = 0.5220 +0.0966, pp = 0.3107 40.0391, p, = 0.7014£0.0573 
and pe = 0.6259 + 0.0603. 


Figure 4.2 Section by the plane (0.01, 0.5, 0.5, -0.5). 
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When sectioned by the plane (a,b,c,d) = (0.5, —0.5, 
0.5, 0.01), our VT gives a picture as shown in Figure 
4.3. Here the points shown are merely the average 
of the coordinates of all vertices in each cell. From 
fourteen simulations of these 50 cells, 104 bonds, 140 
vertices and 187 edges having respectively the coor- 
dinate numbers of 4.1600, 8.2500, 2.6714 and 3.5080, 
we obtain p, = 0.6914 + 0.1424, p, = 0.4533 + 0.1108, 
Py = 0.7760 + 0.0821 and p, = 0.6929 + 0.0897. 


Figure 4.3 Section by the plane (0.5, —0.5, 0.5, 0.01). 


When sectioned by the plane (a,b,c,d) = (—0.7, —0.3, 1, 
0.01), our VT gives a picture as shown in Figure 4.4. Here 
the points shown are merely the average of the coordinates of 
all vertices in each cell. From twenty simulations of these 121 
cells, 305 bonds, 291 vertices and 411 edges having respec- 
tively the coordinate numbers of 5.0413, 9.4426, 2.8247 and 
3.7518, we obtain p. = 0.5413+0.1107, pp = 0.3584 +0.0494, 
py = 0.7077 + 0.0572 and p, = 0.6749 + 0.0470. 


Figure 4.4 Section by the plane (—0.7, —0.3, 1, 0.01). 


For the purpose of doing these simulations, the code for finding percolation has been rewritten 
into a function. To my surprise and delight, the same function works for both the 3-d VT and its 
2-d section. In choosing a sectioning plane it is better if we choose the parameter d small, as the 
plane will then pass close to the origin. Also, choosing a+b+c % 0 seems to make a more wholesome 
section than otherwise. The codes for sectioning work well for oblique planes but do not like planes 
which are parallel to an axis. This shortcoming can be avoided if we make our plane only nearly 
parallel, when we want it to be parallel to an axis. Then to be able to view in a head on fashion such 
planes which have been plotted in three dimensions, we can look from the position (a, b,c), which 
is in effect the vector normal to the plane. The function mentioned above is listed as perc.m and is 
in § A.26. 


§ 4.5 Network percolation 


In the study of networks an important parameter is the coordination number, which is the number 
of neighbours of an element which in the graph theory is usually the vertex. Each vertex or site of 
a graph is connected to each of its neighbouring vertices by a bond, so the coordination number of 
a graph is the number of bonds connected to a vertex. 

Clusters and their various characteristics play an important role in the study of percolation 
of networks. With a material science application in mind, Levy et al (1982) numerically represent 
the shape of a cluster by the shape parameter S, defined as S = b/N or more generally S = 
(1/2) 07_, vi / O31 % where b or i is the number of bonds and N or v; the number of elements in 
a cluster having b or 2 bonds respectively. 


§ 4.6 Percolation statistics in literature 


Rushbrooke and Morgan (1961) studied the Heisenberg and Ising ferromagnetics and found that the 
critical concentration p, is the same for both Heisenberg and Ising problems. Moreover, this value is 
irrelevant to the magnitude of the spin s. By studying point-clusters and link-diagrams, they found 
approximations to the exact values of the percolation probabilities for, the face-centred cubic lattice 
Pe = 0.18, the body-centred cubic lattice p. = 0.22, tthe simple cubic lattice p. = 0.28, and the 
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plane square p, = 0.48. Only vertices of their point-clusters, not link-diagrams, need to have a bond 
between them in order to be neighbours. 

Frisch et al (1961) give the site critical probabilities for various types of lattices in two- and 
three dimensions. In the terminology used here, these are p, of vertices of the networks considered. 
Their results are shown in Table 4.4. 


lattice dimension z pc (N=1,000) p- (N=2,000) 
triangular 2 6 0.487+0.021 0.493 + 0.018 
square 2 4 0.575+0.017 0.581 +0.015 
hexagonal 2 3 0.683 £0.020 0.688 + 0.017 
h.c.p. 3 12 0.204+0.008 7 
f.c.c. 3 12 0.199 + 0.008 - 
simple cubic 3 6 0.325 + 0.023 - 
tetrahedral 3 4 0.434+0.013 0.436+0.012 
ice (quartz) 3 4 0.432+0.013 0.433+0.011 


Table 4.4 Critical probabilities from Frisch et al (1961). 


From their investigations they also realise that the critical probabilities of two homohedral 
tilings may not be the same. For instance, a plane lattice with z = 3 may have the vertex- and edge 
percolation probabilities different from those of the hexagonal lattice, even though the coordination 
number is the same for both. 

Both Monte Carlo technique and series expansion method are means by which one can nu- 
merically study percolation problems. Dean (1963) used Monte Carlo method because, according to 
him, it provides more precise information. In a network of N sites pN of which are occupied, the 
probability of site occupation at the next time step becomes g = p+1/N. The critical probabilities 
he found are shown in Table 4.5. 

lattice size 


12 x 12 24 x 24 48 x 48 
lattice type De +8.d. De +s.d. pe ts.d. pe + 8.d., (lattice size) 
square (s) 0.507 £0.090 0.582£0.032 0.5804£0.018 0.569, (78 x 78) 
triangular (s) 0.435 +0.029 0.494+0.037 0.486+0.017 0.486, (84 x 84) 
square (1% and 2") (s) 0.322+0.047 0.38140.029 0.387+0.014 0.401, (90 x 90) 
honeycomb (s) 0.641+0.061 0.67540.027 0.688 +0.015 0.679, (2 x 84 x 84) 
kagome (s) 0.609 + 0.047 0.643+0.028 0.635+0.020 0.655, (72 x 72) 
four-eight (s) 0.679 +0.039 0.718+0.023 0.732+0.015 0.739, (84 x 84) 
square (b) 0.468+0.049 0.469+0.028 0.492+0.011 0.498, (2 x 96 x 96) 
triangular (b) 0.279+0.038 0.324+0.046 0.329+0.021 0.349 
kagome (b) 0.419+0.058 0.432+0.045 0.449+0.032 0.435, (84 x 84) 
four-eight (b) 0.615+0.050 0.649+0.028 0.675+0.027 0.661, (2 x 78 x 78) 


Table 4.5 Critical probabilities and the standard deviation (s.d.), pe + s.d., for sites (s) and 
bonds (b), from Dean (1963). 


The effect of shape of arrays found by him is such that for array sizes of 6 x 96, 12 x 48, and 
24 x 24 the percolation probabilities are respectively 0.707 + 0.043, 0.594 +0.039, and 0.568 + 0.032. 
The modified second moment of the cluster size distribution for a lattice has been defined as p53 = 
> 0?/(° 0%)”, where o; is the size of the i**-cluster. And the percolation probability for the finite 
lattice has been defined to be the value of p at which Ay3/Ap is maximum, A being an increment 
of one step. 

Tiyapan (1995, KNT3(ii)) finds from 27 simulations on 2-d Voronoi networks of between 70 and 
500 cells the percolation probability of cells p. = 0.507, and from 75 simulations on 2-d Voronoi 
networks of between 100 and 8,600 bonds p. =0.658. For honey comb lattices between 100 and 
3,000 bonds he finds the bond percolation probability p. = 0.640, for Kagome lattices between 300 
and 3,000 bonds p, = 0.517, for square lattices between 200 and 2,000 bonds p, = 0.467, and for 
triangular lattices between 300 and 3,000 bonds p, = 0.341. These results from simulations are 
usually less than the exact values by a few per cent. 

Mecke and Seyfried (2002) find pseudocritical threshold for different types of lattices and then 
extrapolate these by using finite-size scaling laws to obtain the critical probability for the infinite 
network. They use the Hoshen-Kopelman method to find the largest cluster of the system and then 
determine the percolation in fifteen steps by increasing or decreasing p by ae depending on whether 
the system percolates or not. Here Ap = po — pi, pg > p1, such that the system only percolates 
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at po not p;. In two dimensions the network percolates when two opposite sides connect with each 
other whereas in three dimensions the same is true when four opposite sides connect together. The 
pseudocritical percolation threshold is p.(£), where L is the network size. The finite-size scaling 
theorem is p.(L) = p.(oo) + aL’. In two dimensions they find p, =0.59278(4), while in three 
dimensions p, = 0.31162(8). 


In its early days the percolation theory concerned itself much with self-avoiding walks. The 
connective constant or the walk limit is defined as a measure of the connectivity of the lattice as 
In = limyp-+.0(1/n)Inc,, or sometimes as & = In. Shante and Kirkpatrick (1971) gives results 
from both the Monte Carlo and series method as shown in Table 4.6. 


z bb Monte Carlo series method 
Pe Pv Pe Po 

honeycomb 3 1.8484 0.640 0.688, 0.679 0.6527 0.700 
kagome 4 0.435 0.655 0.6527 
square 4 2.6390 0.493,0.489 0.581, 0.569 0.5000 0.590 
triangular 6 4.1515 0.341, 0.349 0.493, 0.486 0.3473 0.5000 
diamond 4 2.878 0.390 0.436 0.388 0.425 
S.C. 6 4.6826 0.254 0.325 0.247 =—0.307 
b.c.c. 8 6.5288 0.178 0.243 
f.c.c. 12 10.0350 0.125 0.199 0.119 0.195 
h.c.p. 12 0.124 0.204 


Table 4.6 Critical probabilities and connective constants, Shante and Kirkpatrick (1971), where 
s.c. means simple cubic, b.c.c. body centred cubic, f.c.c. face centred cubic and h.c.p. hexagonal 
close packing.For critical probabilities, all the values given in four decimal places are exact where 
as those given in three decimal places are numerical. 


§ 4.7 Percolation of n-gons in continuum 


Consider the continuum percolation of regular polygons in two dimensions. Polygons are placed 
on the plane randomly with regards to both their position and orientation. From their position one 
can find the Delaunay triangulation. Using the bonds of the triangulation obtained as reference 
axes, one for each pair of polygons, one then finds the orientation of any two neighbouring polygons 
with reference with the axis connecting them. 


Unlike in the case of Voronoi percolation where the size of each cell depends on the density 
of the generating points, in continuum percolation the size of the polygons has been decided in 
advance. The expanse of the region of interest is then determined relative to this size. The area 
being simulated can then be imagined as lying within an infinite plane where the number density of 
polygons is homogeneous. 


The area of the n-gon is A = n-2- (1/2)yz. 
But y = Rsin@ and z = Rcos6, where R is the y 
radius of the circumscribing circle. Therefore 
we have R = A/(nsin@cos@). The radius of 
the inscribed circle is r = z. 


A list is made of every angle that the rays 
perpendicular to the sides of the polygons make 
with the x-axis in the counterclockwise direc- n 
tion. Then the key procedure of the algorithm is 
to find for each pair the ray which lies closest to 
the reference line joining the two polygons. The an- 
gles that these two rays make with the reference axis 
determines whether or not the two polygons intersect. Figure 4.5 
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OK a The inscribed circle and the circumscribing circle both play a role 

\ = ra y) in deciding whether two polygons form a cluster. Let the radius of the 

ee Sie _ J former be represented by r and the radius of the latter R. If any two 

ail AN May <= polygons are less than 2r apart, then they must overlap. They just 
Sey touch if they are exactly 2r apart. And if the distance between them is 


more than 2R, then they can never touch each other. When the distance lies 
between 2r and 2R there is a probability that they will touch. Whether this 
Figure 4.6 is the case or not depends on the orientation of the two polygons. 


Among all the rays similar to those 
shown in Figure 4.7 one can find the one 
closest to the line between the centres of 
the two polygons, ie 6; and 62 of Figure 
4.8 are a pair of minimum angles for c, 
and cp. Let the centres of these poly- 
gons be c; and cg, and the line joining 
them @jCy. Two such lines obtained 
from the two polygons can fall into on 
of the two cases; they can either be 
on the same side of Cjé), or they are 
on the opposite side of each other as 
shown in Figure 4.8 (a) and (b) re- 
spectively. In either case, the following 
theorem is true. 


Figure 4.7 


Theorem 4.1. The line segments @@, and bby always intersect T{ép. 


Proof. Without losing generality, consider only the case of G@j@». Suppose the above statement is 
false and @j@_ did not intersect C[G2. Then, because @ Gp» is a part of the circumference of a polygon, 
there would always be another ray dy originating from c, which has aja; connected to it, such that 
Ga; intersects €[C2, 1,7 and k being positive intergers. But the two triangles Aca a2 and Ac, aja; 
are identical. Therefore the angle 4; between aja; and @]¢c, would be smaller than @; and we replace 
d, with dy, because the latter is closer to @¢). o 


Furthermore, between 9; and 62 the one which is greater would belong to an aggressive cell, 
while the other one is in a passive state. Suppose that 6; > 02, then c, would be the one which does 
the touching first. Theorem 4.2 states this in a more formal manner. 


Theorem 4.2. If the distance between any two identical regular polygons allows them a certain 
probability p of overlap, such that 0 < p <1, then the polygon whose ray is furthest from the line 
connecting their centres will be the one which causes the two to intersect. 


Proof. Let c, and cz be the two identical polygons in question. As mentioned earlier, there are 
two cases to be considered, which are represented by Figure 4.8 (a) and (b). Because the distance 
between them is more than the diameter of the inscribed circle, when 6; = 62 = 0 the two polygons 
do not touch each other. The two polygons being identical, it is suffice to consider only one of 
them. But because there is an overlap area between the two as shown in Figure 4.7 where there is 
a probability of colliding to occur, and as the distance from a centre c, increases as one goes from 
a midpoint d, of an edge to the next vertex d2 closest to (Co, the greater 6, is, the further away 
from c; is the point of intersection between d,b. and (ep, the more that polygon overlaps into the 
collision hazard zone mentioned above, and therefore the more contribution to the probability of the 
collision. D 


Alternative proof. As shown in Figure 4.8 (c), with 6; > 02 it is the corner a2 which collides with 
the flat of the side b,b2 and penetrates into the second polygon. 
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Imagine the polygons revolving around until they 
collide. At the point of collision if their vertices 
touch, then 6; = 69 and it is a special case where 
both contribute equally to the touch. The case of 
Figure 4.8 (b) where the two angles lie on the oppo- 


site side of each other follows a similar line of rea- 
ain Figure 


soning. 
4.8 The distance between c, and cz is less than the 
diameter of the circumscribing circle but more than 


the diameter of the inscribed circle. 


(c) 


Figure 4.8 When two n-gons practise the Ayudhya 
sword together. 


Figure 4.8 looks as though two n-gons were practising with each other the Ayudhya dab (sword). 

If 6; > 2, then d(@a3) > d (c2b3). Furthermore, a touch or penetration implies that d (¢a3)+ 
d (cob) = d(G@@). If d(G@a3) + d (c2b) < d(Gey) then the two polygons do not touch each other. 
The case when they merely touch at vertices is the case where d (cob) = d (cobg). On the other hand, 
the condition d (a3) + d (c2b3) > d(G@e) does not say anything much because az may just pass 
under bz without touching the side in Figure 4.8 (a), or it may be the case where the two opposing 
sides lean towards the same direction as in Figure 4.8 (b). Therefore, we now have the following 
theorem which is later used in Algorithm 4.6 for finding overlaps. 


Theorem 4.3. Any two polygons p, and p2 touch each other if and only if d(¢ja3) + d (c2b) = 


d (G2). 
Proof. See explanation above together with Figure 4.8. D 
Algorithm 4.6 Find continuum percolation of regular polygons on a plane. 


place polygons randomly; 
find radii of their inscribed circles and the circumscribing circles; 


find coordinates of the vertices; 
find Delaunay triangulation; 
for i = 1 to number of triangles do 


for j = 1 to 3 do 
find distance between the j'"-pair vertices of the i*” triangle; 


if this distance < 2r then polygons overlap do 


record cells overlap; 
elseif distance < 2R then polygons could overlap do 


record cells to find whether overlap; 


endif 


endfor 
find angle of each of the two polygons when viewed from the other 


endfor 
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for every side of every polygon do 

find angles of the normal ray counterclockwise from the positive x-axis; 
endfor 
for every cell which could still overlap do 

find normal rays closest to the lines going towards its neighbours; 


endfor 
complete the overlap check for all remaining cells; 
check percolation; o 
po : 5 
. e Qs 
om OneY . 
a - aa oS 
(O) a ene 
—- _ a bY 
NYE (}) 
1| RO, go 
YOOIY 
SY SD a aneer, oe Nay 
Cf i £ E st ( a hy 
(b) (c) 


Figure 4.9 (a) Randomly placed triangles. Here the program tells us that the largest cluster 
has seven triangles, and the next largest one has six. There are 8 isolated triangles, 5 clusters 
of two, and three clusters of three. According to the program, the two clusters of sizes six 
and two on the upper right corner do not quite touch. (b) The largest cluster in this case 
has got nine pentagons. (c) These 11-gons are also too sparsely placed to percolate. The 
largest cluster in this case has got 8 members. The number of all 11-gons is 40, the same as 
in the two previous cases, and six of which are isolated. There are one cluster of four, four 
clusters of three, and five clusters of two. 

The program in § A.5 uses the theorems above to find whether an aggregate of polygons 
percolate or not. Here each polygon has a unit area. The area of an n-gon is nr? cos(@/2) sin(6/2), 
where r is the radius of its circumscribed circle. By setting this area to one we can find r, which, 
together with the decided orientation of each polygon gives the coordinates of its vertices. The 
program is explained in Algorithm 4.6. 

A percolation threshold normally means the critical percentage by number. But in the case of 
percolation in a continuum where there is no fixed amount of lattice sites to refer to, the per cent 
area covered may be a better candidate for p,. Finding this area can become quite computationally 
intensive it is a problem of finding the union among sets which, for the case of three sets is AUBUC = 
A+B+C-—AB-AC+ABC, for four sets AUBUCUD = A+B+C+D-—AB-—AC—AD-BC 
BD-—CD+ ABC + ABD + ACD + BCD — ABCD, etc. For the general case where there are n 
sets intersecting one another, then, U,, Ai = 0, At — Yona) Aig + Dina) Adie — 1 +n Atom + 
(1D) MOA gay = R(-D“M Ain, where Aja... means Aj;N.AjN A, N---, the subscripts (n, k) is 
the combination "Cx = n!/[k!(n — k)!] and i = ijk--- up to the n* term. 


§ 4.8 Polygon percolation threshold 


Because, unlike lattices however random, the continuum has neither underlying sites nor ver- 
tices and therefore has no definite number of these. To describe the percolation threshold in such 
situation we need to define the some new parameter other than our usual p,., which not only de- 
pends on the number of polygons but also their area relative to that of the domain which they are 
in. We define the threshold area ratio as Definition 4.1. This definition disregards the intersection 
of the polygons, and therefore can only be used to compare networks which have the same type of 
distribution. 

Definition 4.3. For an aggregate of n-gons all the centres of which are Poisson points within a 
space of area A, the threshold area ratio is an = >> ,,|a;|/A, where a; is the area of each n-gon. If 
all n-gons have an equal area, then a, = na/A.: 

The program in § A.5 is transformed into a function which is then reused several times in 
the course of operation of another program, listed in § A.33, which finds a,. Both the program, 
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ppgt.m, and the function, cmpc.m, are here listed concatenated together whereas in practice they 
must separately exist. 


The program first proceeds to find whether percolation occurs for aggregates containing in- 
creasing numbers of pentagons in step of 16. This step size is in fact 2" for some integral k, and 
can be chosen to be optimum for a certain size of the network. At the end of these steps we find 
whether our aggregate percolate. If it does, then we know that somewhere within the last increasing 
step exists the point where percolation sets in. The precise location of this point can then be found 
by doing binary searching. There are altogether k steps of these searches, at the end of each one of 
which there is again a test for percolation. If this test for the last and k*" step is successful, then 
the last picture represents the aggregate which starts to percolate, otherwise one needs to know 
how many consecutive percolation tests have failed up to that point. With the knowledge of that 
number, one can count backward and reach the picture in question. The program ppgt .m fails to do 
this last job, but it has been added, together with the generalisation to 2" steps, to make ppgk.m. 


Figure 4.10 is the sample run on an aggregate of pentagons. Each of the pentagons is one unit 
area, and the space containing them, or rather the space which contains all their centres, is 10 x 10 
in size. 
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Figure 4.10 Finding the threshold area ratio of pantagons. The number of pentagons from (a) 
to (h) increases in step of 24 = 16. In (i), (j), (k) and (1) a binary search proceeds which divides 
the length of the last interval from 16 to 8, 4, 2 and 1 in sequence. This simulation gives a5 = 
1.28. 


The program ppgk.min § A.33 uses k = 4 as explained above. Generalising this into a general 
k gives Algorithm 4.7. Here v;; is the j*" vertex of the i” n-gon. At the completion of the program, 
if it does complete, we have the aggregate containing m4 /2 polygons, and our b' percolation test 
before last has been on this aggregate. 


Algorithm 4.7 Threshold area ratio. 


st 2k; 
Be 2n/n; 
r & [1/(nsin(8/2) cos(3/2))]?; 
while percolation check fails do 
{(x,y)} < find s more centre coordinates; 
{9} < find s more orientational angles; 
for all these new centres and angles do 
vig — (@ + rcos(; + J8),y + rsin(O; + Jf); 
check percolation all the m centres and angles; 
endfor 
endwhile 
(mo,1™m1] < the last step above; 
b+ 0; 
repeat & times 
m4 2 <= (no +1)/25 
check percolation all the first m /2 angles and centres; 
if percolated then 
my — M4123 
b+ 0; 
else 
mo — M4123 
b+ b41; 
endif 


endrepeat o 


Figure 4.11 shows the percolation of n-gons with various values of n, the results of apply- 
ing Algorithm 4.7 above. Polygons up to the twelve-sided dodecagon are shown here. After this 


there are of course the triskaidecagon, tetrakaidecagon, ..., enneakaidecagon, icosagon, icosikai- 
henagon, icosikaidigon, icosikaitrigon, ..., icosiaienneagon, triacontagon, triacontakaihenagon, ..., 
triacontakaienneagon, tetracontagon, ..., pentacontagon, ..., hexacontagon, ..., heptacontagon, ..., 


octacontagon, ..., enneacontagon, ..., hectagon, etc. But all of these closely approaches the circle. 
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(ci) (cj) (cr) 
Figure 4.11 Percolation of n-gons. Triangle, square, pentagon, hexagon, heptagon, octagon, 
enneagon, decagon, hendecagon and dodecagon are represented respectively by t, q, p, h, k, o, 
e, d,n andc. The subscripts i, j and k are respectively the fourth step of the program, the 
over-percolated case and the critically percolated case. a3 = 1.03, a4 = 1.17, as = 0.89, ag = 
0.99, a7 = 1.28, ag = 1.038, ag = 1.28, aig = 1.07, a1, = 0.91 and aig = 0.95. 


§ 4.9 2-homohedral tilings 


There are altogether 39 types of 2-homohedral tilings, ie. those which have vertices of the same 
valence. Of these, there are 26 types of valence 3, 10 of valence 4, and 3 of valence 5 (cf Griinbaum 
et al, 1987). The code and the data used in the simulations are in § A.6. In the program, the variable 
o contains the vertex numbers which must be ordered by scanning the basic tile from left to right 
and gradually from bottom up. Unless the variable follows this ordering strange results will occur. 
The variables m and n are respectively the z- and y-coordinates of the vertices. 
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Figure 4.12 The thirty-nine 2-homohedral tilings and their covering lattices. 


Once while doing the above simulations I came across an interesting batch of the 44[34]72[3"|n 
2-homohedral tiling where the percolations forward and backward in time gave the same value of 
pc, that is there is an approximate time symmetry. In other words the p,’s of both the first and 
the second population are symmetric to each other. Each percolates when its population reaches 
756 out of the combined total of 1070, which identically gives p. = 0.7065. Large non-percolating 
clusters quickly disappears after the onset of percolation as shown in Figure 4.13. In Figure 4.13 
(b) the size of the second largest cluster is reduced to one between the population II of 941 and 961, 
which correspond to the density of 0.8794 and 0.8981 respectively. Interestingly from the population 
of 961 until 1070 there is only one cluster despite the fact that there are still more than one hundred 
population I. Here population IT always refers to any population which percolates or which is being 
observed with regard to percolation. On the other hand in Figure 4.13 (b) the second biggest cluster 
is the smallest cluster from population 898 to 1037. Again from the population of 1038 until the 
maximum 1070 no second clusters appear. 


10° r r r T r 10° 
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Population II cluster size 
Population II cluster size 


L L L L L L L L 
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10 600 
Population Il number 
(a) (b) 
Figure 4.13 The biggest, the second biggest, and the smallest cluster from a batch of the 
44[34]72[3"]11 2-homohedral tiling where the percolation happens to be symmetry in time, that is 
(a) and (b) are time symmetric to each other. The vertical dotted line is where the percolation 
point occurs in each case. The the heaviest to the lightest lines are respectively the progression 
of the clusters which are biggest, second biggest, and smallest at each population density. 


The percolation considered is that which percolates from the lower bound to the upper bound, 
both of which are fixed. Because the networks must be constructed such that there would be 
duplicates of neither vertices nor links, there must be cells some vertices and edges of which belong 
to one or more of the basic tiles adjacent to them. As a result, those basic tiles on the boundary 


0 600 
Population Il number 
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sometimes have a few of their cells missing, and depending on whether the basic tile in question is 
very complex or simple the position of the missing cell may merely lie on its boundary or may lie 
nearly halfway into its interior as is the case with the 43[34]83[3°]; and the 43[34]83[3°]m_ tilings. 
The upper and the lower boundaries, which decide the percolation, are chosen as straight lines 
parallel to the vertical axis and lying in the direction towards the interior of the network and at a 
fixed distance from the the maximum and the minimum horizontal positions respectively of all the 
vertices in the network. This fixed distance is defined for each network to be a fixed ratio of its 
overall width. The value of 0.05 had been used which was later changed to 0.1 due to a problem 
which occured while simulating the 42[3*]82[3°] tessellation, missing cells formed a continuous band 
disconnecting all the cells on the boundary from the rest of the network, the upper boundary lay 
within this band when shifted by five percent from the upper limit but not when shifted by ten 
percent. 


The tilings in Figure 4.13 have the basic units shown in Figure 4.14. 
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3a[4°]53[45] 31[4°]62[4°]1 31[4°]62[4°]n 32[4°]64[4°] 33[4°]66[4°] 


31 [49]84[43] 31 [5°]40[54]y 31[5%]40[54] 3o[5°]44[54] 


Figure 4.14 The unit cells of the 2-homohedral tilings. 


In Table 4.7 there are six types of percolation probabilities. Two cells can be neighbours to 
each other under the criterion of having either one or two vertices in common. This gives rise to two 
different values of p,’s for cells, and another two for bonds. These can be termed 1-neighbours and 
2-neighbours respectively. In Physics and other physical sciences 2-neighbours are more important 
because the greater the contact area is the stronger the bonding. This is not always the case in other 
fields, for example in Sociology a smaller point of contact sometimes implies the less duplications 
of neighbours’ neighbours, and thus the larger the resultant networks. Because there is a cost for 
maintaining the connections, the smaller each contact point is and the larger the network the better, 
therefore here 1-neighbours seem to be more important. 


1. 33[37]71[37] 


Cell Bond cell bond vertex edge 
0.471420.0667 0.3041=50.0272 0.7212-20.0444 0.6771=20.0381 
1.1551 x 10-6 —3.3740x 10-6 —1.1124x10—4 —5.2007x10—5 
4.0885 x 10-5 1.1316x 10-6 1.9314x10—5 9.5566x 10-6 


5/99] 5 4195], 4[256], 4[470],] [4[349], 4[657], 
5.2525 5.4667 10.4149 2.7266 2.7957 3.7077 3.7869 
12[391] 12[894] 12[1275 

5.6215 2.8523 3.8494 
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2. 33[3°]93[3°]n1 


Cell 


0.4966=00.0655 


— 8.2288 


x10—5 


3.7765x10—5 


Cell 


6[412] 
5.5243 


0.4695-=50.0858 
—3.5574x10—5 
1.4884x 10-4 


Cell 


0.4903-L0.1015 
1.2205 x 10-4 
2.1262x 10-4 


Cell 
0.46074 


C0.0586 


7.9409 x 10-5 
3.3840 x 10-5 


Cell 
0.52862 


C0.0511 


5.4350 x 1079 
1.7402x10—5 


Cell 
0.52062 


C0.0551 


—9.0646x 10-6 
1.7837x 10-5 


Cell 
0.47264 


C0.0684 


1.0814x10—4 
8.0038 x10—5 


Cell 


0.49582 
1.3909 
8.0344 


-0.0796 
x1o—-4 
x10-5 


Bond 


0.2734=-0.0206 
—5.5513x 10-6 


cell 
0.504244 


—1.0585x10—4 


bond 


0.3119=-0.0417 
4.8211x10—5 


C0.0574 


5.9281 x10—7 2.1108x10—5 7.5982x 10-6 
6[831],  6[1138] |[2[305], 6/412] | [2[771],  6[1054] 
12.0144 12.1564 5.0557 5.1165 11.0298 11.1556 
4 8 
3. 44134]8,[3°| 
Bond cell bond 
0.3044=-0.0429 
2.8967x10—5 
6.8367 x 10-6 
4[748],  12[888] 
10.5294 10.5946 
3 8 
4. 33[3°]82[3°] 
Bond cell bond 
0.292600.0281 
—1.8541x10—5 
1.9682x 10-6 
4[812],  12[1313] 
11.0788 11.2749 
3 9 
5. 33[8°]93[3" Jin 
Bond cell bond 
0.24832-0.0258 
4.8742x 10-6 
7.4553x10—7 
4[957|,  12[1188] 
12.0878 12.1785, 
4 10 
6. 43[3°]106[3" 1 
Bond cell bond 
0.2991-0.0319 
—8.0790 x 1076 
2.7687 x 10-6 
4/803] , 
11.5691 
9 
7. 33[39]93[3°]1 
Bond cell bond 
0.2987=-0.0308 
7.4580x 10-6 
1.7780x 10-6 
4/699],  4[1313], 
11.9371 12.2224 
6[2119] 
12.3870 
8. 45[34]104[3!°] 
Bond cell bond 
0.2905=-0.0633 
5.5163x10—4 
1.3835x 10-4 
4[1537], 
11.6695 
4 10 
9. 43[3 ]10¢[3 lu 
Bond cell bond 


0.3020=20.0381 
3.9193x10—5 
6.3920x 10-6 

4[866] , 

11.5727 


vertex 


0.7628=-0.0228 
1.0042x 10-6 
7.3299x10—7 


4/415), 
2.8193 


2[692], 
2.8584 


6[919] 


2.8770 


vertex 


0.7194=-0.0209 
—5.8481x10—6 
6.5739 x10—7 


edge 


0.7094 -10.0435 


—2.759' 


ox10—5 


7.7345x 10-6 


edge 


0.6407=20.0296 
2.8990x 10-6 
2.0123x 10-6 


vertex 


0.74470.0329 
2.1925 x10—-7 
2.8427 x10—6 


edge 


0.7019=£0.0270 
—6.4679 x 10-6 
8.9933x 10-7 


vertex 


0.7482 


-0.0264 


2.4849 x10—6 
1.2161 x10—6 


edge 
0.70492 


£0.0195 


—1.0461x10—5 
1.0354x 10-6 


vertex 


0.77892 


C0.0386 


8.7871x10—6 
4.7653 x10—6 


edge 
0.71734 


C0.0391 


—1.2114x1074 
1.8197x 10-5 


vertex 


0.75322 


C0.0296 


—8.9673 x10—6 
1.3152x10—6 


edge 
0.70634 


C0.0222 


—4.5643 x10—6 
7.3255x 10-7 


vertex 


0.73412 


C0.0316 


—8.0603 x10—6 
2.0508 x10—§ 


edge 
0.69574 


£0.0303 


7.5952x 10-7 
1.6121x10—6 


vertex 


0.74902 
1.0834 
1.4492 


-0.0461 
x10-4 
x10-5 


edge 
0.67514 


C0.0199 


6.2752x 10-7 
3.7299x10—7 
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10. 45[34]82[3°] 


Cell Bond cell bond vertex edge 
0.5102=50.0568 0.3183=20.0405 0.7202=-0.0493 0.6953=-0.0307 
4.1767x10—5 —9.0691x10—6 —1.1640x10—4 1.5917x10—5 
3.0853 x 10-5 7.8890x 10-6 2.1475x10—-5 1.8233 x 10-6 


411040], 
10.5096 


11. 43[34]83[3°]; 


Cell Bond cell bond vertex edge 
0.4577=0.0761 0.29102-0.0421 0.6887=-0.0279 0.6532=-0.0340 
—2.5782x10—5 —1.6419x10—5 8.0593x10—7 — 2.5506 x 10-5 
8.4655 x 10-5 7.8109x10—6 7.8808 X10—7 3.2110x10—6 


4978] , 
10.5562 


12. 43[34]83[39}n1 


Cell Bond cell bond vertex edge 
0.5030=20.0603 0.316620.0334 0.7385=20.0255 0.6718=-0.0369 
2.6397x10—5 —1.3936x10—5 —2.6605 x10—6 1.8952x 10-6 
4.0377x10—5 2.7675x 10-6 9.6100x10—7 3.1533x 10-6 
4[219], 4[419],][4[582],  4[1152], 
5.3151 5.4988 | |10.2818 10.5799 


13. 43[34}83[3° Jin 


Cell Bond cell bond vertex edge 
0.4931=20.0560 0.2927=20.0320 0.7294-20.0348 0.6730=20.0289 
1.5620 x 10-4 3.0191x10—5 4.5084 x10—6 —2.7067x10—5 
3.1433 x10—5 6.4721 x 10-6 2.5226 x10—6 3.4615x 10-6 


41133}, 
10.6019 


14. 44[34]72[3 Jn 


Cell Bond cell bond vertex edge 
0.4761=20.0432 0.317420.0396 0.7144=20.0228 0.644020.0371 
—7.3074x10—5 3.7498x 10-5 7.2116x10—6 2.2191x10-5 
9.5057 x 10-6 6.2114x 10-6 6.1119x10-7 3.8797x 10-6 
8[345], 12[486] |/8[961],  12[1370] 4[767], 12[1063] | |4[1070], 12/1500] 
5.5710 5.6379 10.0208 10.1241 2.7901 2.8222 3.8579 3.8800 


15. 33[3°]12,[3"7] 


Cell Bond cell bond vertex edge 
0.5077=C0.0896 0.2640=-0.0333 0.8273=-0.0209 0.7269=-0.0313 
4.5180 x10—4 2.5513x 1076 1.3947 x10—5 —2.2183x 1076 
3.2202 x 10-4 3.2552x 10-6 1.6592 x 10-6 1.6715x 10-6 


4529] , 
14.2949 


12[1343] 


16. 44[34]72[37]1 


Cell Bond cell bond vertex edge 
0.4774L0.0624 0.31400.0390 0.730520.0263 0.650420.0349 
—1.6958x 1074 3.5112x 1075 7.4362 x 107 § 1.7310x 1075 
4.8309 x10—5 9.6223x10—6 1.0424 x10—6 3.0145x 10-6 
8[403], 12/554] | [8[1100], 12[1534] 471050],  6[1421]] [4[1296],  6[1 764] 
5.4591 5.5379 10.0236 10.1213 2.4686 2.4828 | [3.8704 3.8889 


17. 54[3°]74[37]1 


Cell Bond cell bond vertex edge 
0.5112=50.0412 0.3406=-0.0207 0.7029=-0.0507 0.6645=-0.0300 
—5.3542x10—5 2.1641x 10-6 —5.7814x10—5 9.3355x 10-6 


8.4833 x 10-6 5.0113x10—7 1.4169x10—-5 1.5205x 10-6 
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18. 3,[4°]51[4°] 


Cell Bond cell bond vertex edge 
0.4003=20.0534 0.2368-20.0198 0.578620.0400 0.4595-50.0398 0.6184 50.0356 0.5109=20.0531 
—1.4626 x10—5 —5.2284x 10-6 3.9800x10—5 —2.9462x10—5 —1.2175x10—6 2.0112x10—5 
1.5617x 10-5 3.5502x10—7 5.2273x 10-6 5.1123x 10-6 2.6728x 10-6 1.7919x 10-5 
4[508] , 6[795] | [4[1872],  6[2979] ] |4[508], 6[795] | [4961], 6[1521] | |4[569], 6[871] | [4[1024], — 6[1600] 
7.37kk01 7.4943 | [14.1731 14.3411 | 13.7835 3.8264 6.1748 6.2406 3.5993 3.6739 5.6738 5.7387 
19. 3)[49]51 [4°] 
Cell Bond cell bond vertex edge 
0.4065-20.0646 0.2309=20.0241 0.5898-L0.0633 0.439420.0536 0.642320.0381 0.5049=20.0271 
— 9.8055 x 1075 2.2560x 10-6 —1.8568x10—6 —1.3300x 107-4 —5.1659x10—5 —9.8556 x 10-7 
3.7012x10—5 6.6307x10—7 3.3297x10—5 3.6224x10—5 7.7919x10-6 8.8048x10—7 
4[512], 6[648] | 4[1954],  6[2486]][4[512), 6[648] | [4992],  6[1260] ||4[545], 6[685] | [4[1024], —_6[1296] 
7.6328 7.6728 14.2027 14.2912 | 13.8750 3.8889 6.2016 6.2349 3.7578 3.7839 5.8145 5.8349 
5 8 
20. 53[3°]8¢6[3°}i1 
Cell Bond cell bond vertex edge 
0.4631=20.0471 0.31642-0.0255 0.6800=-0.0384 0.6587=-0.0246 
x10—6 —9.9185x10—6 1.5506x 10-6 


4.4523 x10—6 6.3039x 10-7 


x10—-6 
4[1906], 
10.1815 


21. 53[3°]86[3°]m 


Cell Bond cell bond vertex edge 
0.5052=20.0474 0.32112-0.0274 0.7145=-0.0299 0.6469=-0.0242 
—1.6178x 10-6 —9.4980 x 10-6 8.1493 x10—6 —1.0841x10—5 
1.0710x10—-5 2.4172x10—6 1.7096 x10—6 7.2095x10—7 
8[351], 12[814] ]{8[971], 12/2316] 4[804],  6[1782] | [4[ii44],  6[2580] 
5.5328 5.6904 10.0062 10.2366 2.8458 2.8956 3.8531 3.9023 
22. 52[3°]1212[3!2] 

Cell Bond cell bond vertex edge 
0.489720.0718 0.3037=20.0257 0.7384=20.0258 0.6694-20.0297 
—1.4823x10—-4 —1.3145x10-5 3.3765 x10—6 1.5898x10—5 
1.0656 x 10-4 1.4342x 10-6 8.4513x10—7 2.4608x 10-6 


23. 54[3°]74[3 Ja 


Cell Bond cell bond vertex edge 
0.5001=50.0342 0.3327=-0.0329 0.6973=20.0235 0.6688=L0.0303 
1.7585 x 10-5 4.4678x10—5 —7.2307x10—6 2.4051x10—5 


3.4856 x 10-6 5.2573x 10-6 5.9288 x10—7 2.0921x 10-6 


8[371], 12/541] | [81027], 12[1519] 4[835],  6[1194] |[4[1195],  6[1722] 
5.5364 5.6155 9.7254 9.8328 2.8623 2.8844 3.8661 3.8885 


24. 53[3°]86[3°]1 


Cell Bond cell bond vertex edge 
0.4867=00.0459 0.3265-0.0237 0.7108=-0.0200 0.6720=-0.0294 
—1.0891x 10-4 —3.7362x10—6 2.8929 x10—6 6.2044x 10-6 
3.3218 x10—5 7.9808x 10-7 4.2792x10—7 3.1424x 10-6 
6[214], 8[535],] 4[577],  8[1502], 2[506],  4[1192],] /2[731], —-4[1 762], 
5.3925 5.6150 | |9.8024 10.1278 2.8261 2.8859 | 13.9398 4.0204 
12[790] 12[2245 6[1728] 6[2573] 
5.6835 10.2254 2.9051 4.0459 


25. 45[34]12,[31?| 


Cell Bond cell bond vertex edge 
0.5117=-0.0960 0.2637=-0.0295 0.7861=-0.0215 0.727320.0227 
1.0580 x 10-4 —4.6675 x 10-6 7.1685 x10—6 1.1750x 10-6 


2.1578x10—4 1.9508x 10-6 4.3496 x10—7 7.1142x10-7 


26. 45[34]1812[319] 


Cell Bond cell bond vertex edge 
0.5249-50.1168 0.24932-0.0632 0.8497=L0.0234 0.7743-0.0483 
1.1603 x 10-3 7.1820x 10-4 5.2173 x10—6 8.6384x 10-5 


6.1479x 10-4 2.0594x 10-4 6.1183 x10—7 1.2578x10—5 


4]1101], 
15.8656 
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27. 53[3°]73[3"]r 


Cell 


0.4882=20.0494 
—4.4634x10—5 
1.5462x10—5 


Bond 


0.3335-20.0310 
—1.2425x10—5 
2.4243x 10-6 


4[170], 4/399] | [4/455], 4[1112] 
5.3529 5.5739 9.4110 9.7392 
Cell Bond 


0.4345-50.0489 
8.3482x10—5 
4.0009x 10-5 


Cell 


0.4099-=50.0264 
4.7650x 10-8 
1.0783x 10-6 


Cell 


0.432250.0512 
5.8373x10—5 
1.9444x 10-5 


Cell 


0.5029-50.0759 
8.6335x10—6 
1.4737x 1074 


Cell 


0.4339=50.0259 
4.6901x10—7 
6.9917x10—7 


2161] 


0.2930=-0.0243 
—9.3058x 10-6 
1.1454x 10-6 
6[1670] 


9.8137 11.3257 


Bond 


0.2271=-0.0154 
2.1524x10—6 
1.2777x 10-7 

2[1407],  6[1863] 

14.0640 14.1857 


Bond 


0.2344-50.0171 
—1.5204x10—6 
2.2157x10—7 


21831], 6/2927] 
14.1038 14.2890 
Bond 


0.2504 -£0.0233 
—3.9666x 10-6 
8.6134x10—7 


cell 


bond 


28. 53[3°]73[3 In 


cell 


0.4785=20.0354 
4.7566x 10-5 
7.5326x 10-6 


2[167], 


2[61], 


29. 30/4 


cell 


0.5687-20.0469 
2.6238 x 1075 
8.5056x 10-6 


2[277], 21379], 


bond 


0.3239-20.0187 
—4.1986x10—6 
3.8717x10—7 


315 9[4°]1 


bond 


0.5028=-0.0283 
1.1978x 1075 
1.1260x 10-6 


30. 32[47]52[4°]n 
bond 


cell 


0.5638=20.0590 
4.8226x 10-6 
1.5253x 10-5 


2[276], 


21496], 


cell 


0.623820.0601 
—1.3226x 10-4 
3.8526x10—5 


0.4749-20.0440 
—5.1510x10—5 
1.0399x 10-5 


3)53[4°] 
bond 


0.5029=-0.0359 
1.1319x 10-4 
1.9717x 10-5 


vertex 


0.7149=20.0273 
1.6306 x 10-6 
9.8464x10—7 


edge 


0.6495-£0.0278 
5.9525x 10-6 
1.4733x 10-6 


2[413], 2907] | [2[574], _2[1293] 
2.7797 2.8512 3.8084 3.8716 
vertex edge 


0.6900 20.0283 
—7.8267x 10-6 
1.5860x 10-6 


vertex 


0.6228=-0.0346 
—1.1730x10—5 
5.5785x 10-6 


vertex 


0.6260 20.0368 
—2.9282x10—5 
6.7165x 10-6 


vertex 


0.6645-00.0414 
1.1038x 10-4 
1.6270x10—5 


0.64212-0.0289 
3.4481x 10-6 
1.8973x 10-6 


edge 


0.5022=-0.0343 
3.4916x10—5 
3.2050x 10-6 


edge 


0.5468-20.0270 
8.1058x10—6 
1.6103x 10-6 


edge 


0.5801=£0.0306 
—3.3375x10—5 
3.6510x10—6 


2[1139], 
13.8191 


Bond 


6[1743] 
14.0425 


0.23352 


£0.0200 


6[161], 
3.5404 


6[473] 
294 


cell 


32. 31[4°]62[4%] 


0.54682 


C0.0321 


bond 


0.42352 


C0.0152 


vertex 


0.64182 


C0.0638 


edge 


0.57152 


£0.0431 


—3.2315x 107-6 
3.6617x10—7 


3.4990x 1075 
4.5289x 10-6 


2.5190x 1076 
1.5119x10—7 


2760],  6[1190) ||2[2860], 6/4534] ||2[760], 6[1190)] | [21457],  6[2301] 
7.5263 7.6202 15.2503 15.4027 3.8342 3.8672 6.7138 6.7718 
33. 31[4°]62[49]n 
Cell Bond cell bond 


0.4740=50.0681 
—1.5376x10—5 
4.0570x 10-5 


Cell 


0.454020.0687 
7.8079x10—5 
4.6624x10—-5 


0.2803=20.0276 
—1.1319x10—5 
1.3703x 10-6 


Q[1742], 6/2656] 
14.8576 15.0715 
Bond 


0.2265=-0.0369 
—1.5535x 1075 
2.8752x 10-6 


0.6333=200.0534 
—1.4476x10—4 
2.7149x 10-5 


0.4785=20.0260 
1.8163x 10-5 
1.3662x 10-6 


—1.1180x 1074 
4.6849x 10-5 


2[833], 6 [1281] 
3.6879 3.7471 


vertex 


0.6961 =20.0552 
2.5890x10—5 
2.2393x10—-5 


2[532], 


2[284], 


4.9165x 1075 
7.7330x10—6 


2[1536],  6[2400] 
5.7513 5.8008 


edge 


0.6132=20.0465 
6.0868x10—5 
1.1922x10—-5 


34. 32[47]64[4°] 


cell 


0.6160=20.0685 
—1.8755x10—4 
5.5163x10—5 


bond 


0.4908=20.0346 
—1.2315x10—5 
2.5752x 10-6 


vertex 


0.6320-£0.0532 
1.0557x 10-4 
2.0627 x10—5 


edge 


0.5527-20.0434 
1.0264 x 10-5 
1.0088 x 10-5 


2[1605], 
14.9533 


6[2026] 
15.0642 
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35. 33[4°]66[49] 


Cell 


0.4042--0.0580 
3.6150x10—5 
1.8890x 10-5 


Cell 


0.4525-50.0780 
7.1913x10—4 
1.4678x 10-4 


Bond 


0.2281 =-0.0277 
2.6810x10—6 
7.5233x10—7 


2]1782],  6[2418] 
14.8911 15.0488 
Bond 


0.2122-0.0230 
—9.6553x 1076 
8.2295x10—7 


cell 


0.5958=00.0567 
—4.2796x10—5 
1.4656x 10-5 


cell 


0.5237=20.0627 
—1.0143x1075 
3.2798x10—5 


bond 


0.4698-20.0313 
3.1642x10—5 
4.6775x 10-6 


bond 


0.3807=-0.0401 
3.5308 x 10—6 
5.4378x 10-6 


vertex 


0.6353=-0.0449 
—2.0589x 10-5 
8.7591x 10-6 


21543], 


2[267], 


vertex 


0.7010=-0.0470 
—7.5034x 1075 
1.2495x 10-5 


edge 


0.5210-£0.0370 
7.9238x 10-6 
2.4104x 10-6 


edge 


0.6086=-0.0679 
1.6014x10—4 
3.8208x10—5 


[1794], 
16.9197 


6[2632] 
17.1071 


[288], 


2[561], 


Cell 


0.4047=-0.0637 
6.1235x10—5 
3.1306x10—5 


Cell 


0.3457=00.0484 
8.5340x10—5 
2.5074x10—5 


Bond 


0.2185=-0.0212 
5.5289x 10-6 
5.4838x 10-7 

2[2374],  6[3273] 

16.0126 16.1577 


Bond 


0.1820=20.0200 
—5.3816x 10-6 
6.0154x 10-7 


37. 3:[5°]40154] 


cell 


0.6180=20.0491 
6.5660x10—5 
1.1080x 10-5 


cell 


0.6226=20.0443 
—5.2469x 10-5 
1.2122x10-5 
[213] , 2[380], 
3.1632 


bond 


0.4930-20.0327 
1.4895x 10-5 
2.1733x10—6 


bond 


0.5695-20.0477 
—1.2139x10—4 
1.8902x 10-5 


vertex 


0.5516=-0.0469 
4.3305x10—5 
1.3841x 10-5 


2[281], 2[613], 


vertex 


0.5527-20.0393 
—2.9005x10—5 
1.0732x 10-5 


8169] , 2[289], 


edge 


0.4101-0.0301 
1.2816x 10-5 
2.1196x 10-6 


edge 


0.4162-0.0419 
—3.2254x10—5 
7.5538 x 10-6 


Cell 


0.3672=00.0405 
— 6.6453 x10—5 
1.1174x10-5 


Bond 


0.1766=£0.0243 
—3.2608x10—6 
1.0035x 10-6 


cell 


0.6265£0.0398 
7.1192x10—5 
7.0957 x 10-6 


bond 


0.5822-F0.0505 
2.3922x10—5 
2.8516x10—5 


vertex 


0.5413=200.0515 
—9.8715x10—5 
1.9319x 10-5 


edge 


0.4101-0.0536 
9.3126x10—5 
3.1517x10—5 


[732] ; 


2[1 212], 


4[1 68], 


2[270], 


[141], 


4.216], 


from top to bottom, mean + standard deviation, second moment, and third moment. Similarly 
each of the items in each box contains the number of runs [the number of network components] 
and the coordination number. 


The advantage of sparse- over dense networks described by Burt (1992) is an idea similar to 
that of the benefit of decentralisation and local autonomy in politics. According to him, size still 
matters but the cost of maintaining the size is also important. The idea can be very general; dense 
networks are virtually worthless monitoring devices, while sparse networks give more information 
benefits. Taking opportunity costs into account, the rate of return of a dense network is lower than 
that of a sparse network. The idea can be applied to a strategic network expansion which is crucial 
for sales persons and jobs hunting alike. 

Something similar to this comes up in a variety of fields. In Sociology it is sometimes called 
the strength of weak ties. A tie between two things is weaker the less each of them has in common 
with the other. It is well known, i.e. the origin unknown, that best friends are usually those who 
have least in common with each other. It would be interesting to trace this idea back to where it 
first appeared in literature. 

The study above shows that the percolation probability is not a function of the coordination 
number of the network alone. I suspect that it is a function of both this and the coordination number 
of the dual lattice of the network, that is p. = f{(x,,2,). In other words, p, may be a function of not 
only the connectivity but also the tortuosity of the network. 
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§ 4.10 Cosmology 


The universe is made up of superclusters of galaxies, each of which contains local clusters. Each 
local cluster contains galaxies, each of which contains stars. When stars explode, gas bubbles are 
formed which expand to meet one another, forming walls of materials in the manner of the Voronoi 
tessellation in three dimensions. Within the planets around each stars there exist yet other structures 
which can be similarly represented by the Voronoi tessellation. 

Superclusters are like endless cobwebs covering the space. The Andromeda galaxy, the Sculptor 
group, the Virgo cluster, and the M81 group are examples of the members of the Local Supercluster 
in which our Milky Way galaxy reside. The Pisces-Perseus supercluster is 250 million, and the 
Hercules supercluster 500 million light-years away from us. 

Assuming the big bang origin of the universe, the defects originated in the phase transition 
at the Big Bang could play an important role in the formation of galaxies and their clusters and 
superclusters. Tom Kibble predicted this in 1976 (cf Croswell, 1995), and gave three types of 
possible defects, namely the point-like monopoles, the line-like cosmic strings, and the plane-like 
domain walls. 

Van de Weygaert (1989) studied the statistics of 3-d Voronoi tessellations for the purpose of 
understanding the structure of the universe. Here neighbours are defined by a common polygonal 
face whereas full neighbours have an additional requirement that the line joining nuclei intersects 
this common face. Taking a linear section through d-d, d > 2, the mean length is 


ar (4-3) (44) 


(A) = ETON SCC see 
(d—1)lp2ar (4 +1) \*71T (2-4) 


, (20)iy 


where p is the number density, in other words nucleus density. 

In three dimensions we need at least two points in order to establish the distance of an object 
by a geometrical method. Because the cosmic scale is so large compared with the scale of the solar 
system and Earth, all the observation points we may choose become the same and only one point 
in practice. Therefore the determination of the distance to extragalactic objects, which is crucial 
in cosmology, can only be achieved by the various means of observations from a single point in 
space. This is a limitation that has caused the establishment of distance to be a hot, long-standing 
controversial issue. As the universe is expanding with the receding velocity of objects increases with 
their distance away from us, and as this velocity can be accurately determined by the amount of 
redshift in the light of these objects observed, the study of distance becomes the study of a single, 
universal constant called the Hubble constant. This is only a constant in theory not in practice. All 
numerical values of the Hubble constant used in literature are fictitious to some degree. If its exact 
value is known, then the accurate distance to any cosmic object is simply its radial velocity, i.e. 
its velocity away from us, divided by this Hubble constant. The reciprocal of the Hubble constant 
is called the Hubble time. The importance of the Hubble constant, together with the difficulty in 
finding it, make the observational and the theoretical parts of cosmology inseparable from each other. 
They are also the cause of the proliferation of the modern literature on Astrophysics. One example 
of studies of the extragalactic distance scale is that reported by de Vancouleurs (1993) which includes 
among the objects studied the following objects belonging to the Local Group, namely the Large 
Magellanic Cloud, M31, M33, NGC 3109, as well as the objects beyond the Local Group, namely 
NGC300, 1C4182, NGC2403, NGC5128, M81, M101, M104 (NGC4594), and NGC4571. 


§ 4.11 CCTV, forest fire, the navy and porcupines 


Voronoi graph can be used to aid the decision on the locations of Closed-Circuit Television (CCTV) 
cameras. CCTVs have proved to be very effective in reducing the number of crimes in Manchester 
as well as other places within the UK. In order to be effective, it is important to have no large gaps 
within a control area. Strategically placed, the cameras will give the best cost-effective solution 
while covering the largest possible area. 

Within crucial areas cameras could be positioned in such a way that their ranges of operation 
intersect without gaps. The points of intersection represent the vertices of the Voronoi regions centred 
around two or three neighbouring cameras. The design objective is to produce the maximum covered 
area using the least number of cameras. 

Outwards from each crucial area extends space not being covered by cameras. Similarly around 
other crucial areas surrounding this area also extend such unobserved space. The next step in the 
design is perhaps to make sure that people walking from one covered area to another are safe. To 
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do this one can consider the whole covered areas as centres of yet another set of Voronoi tessellation 
orders of unit larger in size. The purpose now is to make sure to position cameras at the positions 
of vertices of the Voronoi regions of this second set of tessellation. 


Both the first and the second tessellations will be subjected to geographical constraints. The 
former would be predominantly influenced by the shape and the position of buildings and sur- 
rounding structures, the latter by the paths connecting the covered regions. Here a vertex of these 
geographically-constrainted Voronoi regions could be a point which is equally far from three covered 
regions. Another set of positions to be considered for placing CCTVs is one of those points equally 
far from two covered regions along the shortest path between them. 


At a graduate development programme on 16” November 2001 I gave a presentation for 
ten minutes on a title division of space. In it I described the current project of mine as concerning 
porous media and division of space in general. The part presented is from a section for miscellaneous 
applications. Mathematical models of partitionings, for example the Voronoi tessellation pioneered 
by G. F. Voronoi and G. L. Dirichlet in the 19" century, help towards understanding physical 
Euclidean and non-Euclidean world. Procedures based on these models help solve many physical 
and strategical problems. When I was fourteen years of age I was trained in a course for those who 
volunteer to fight forest fires in the area around the Sudeb mountain, close to where I used to live. 
Although I never used it afterwards, the knowledge obtained from that course has given me the 
interest in forest fire fighting. 


To mention but one of the relevant topics contained in that 
course, one can imagine a forest fire station as being a nucleus 
of a Voronoi partition of area under the protection of that 
stations. Here the distance from a station to any point on the 
forest ground is not a straight line but the time it takes a fire 
fighter from that station to reach the point in question. This 
time is affected by the tortuosity of the path, in other words 
the degree of winding, as well as the difficulty of the climb. 
There normally is an observation tower at each fire station, 
so we may safely assume that every nucleus is equipped with 
one. Then the exact location of a fire observed can be located 
on the map from the intersection of lines of sight drawn from 
the neighbouring towers towards the direction of the fire. 


Figure 4.15. Range of observatory towers. 
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These lines of sight only give the position of 
the fire on the map. The distance measured 
along such a line does not necessarily corre- 
spond with the Voronoi distance which has a 
unit of time. The area covered by a tower de- 
pends on its height; a higher tower can report 
a fire occurring at a further distance. The 
area of covering is circular if the forest is a 
flat plain. But such a case is rare and in gen- 
eral this area is a distorted circle the degree 
of distortion of which depends directly on the 
contour of the land around it. In Figure 4.15 if 
his the height of a tower and r its radius, then 
rT, <2 implies hg > hy. Suppose the dashed 
lines are boundaries of Voronoi regions, solid 
boundaries the observation range of the tow- 
| ers, and straight lines from nuclei the lines of 
| sight. In Figure 4.16 the Voronoi boundaries 
ee are not straight lines as a result of tortuosity. 
Observations from only two towers suffice for 
finding location of the fire, but an observa- 
tion from the third tower would reduce the 

probability of an error. 


Figure 4.16. Tortuous boundaries. 


T also talked about strategic locations. At a job presentation by Tesco we learnt that the com- 
pany uses the Geographical Information System for placing their new stores in strategical locations. 


The planning department at Tesco uses 
the so called the porcupine diagram which 
looks like the rays of light emanating from 
/ a star or the quills from a porcupine. The 
rays radiate from a supermarket and their 
length represents the drive times required 
A to get to the store from the other end of 
the line. Factors affecting the locating of 
new stores include for example drive time, 
competition, logistic, parking, and traffic. 
road The porcupine diagram superimposed on 
<< Spatial a map. Triangles are Tesco supermarket 
oS pas stores while circles are those belonging to 
e) their competitors. The spatial boundary 
of the Voronoi diagram is in general dif 
ferent from its drive time boundary, as 
shown in Figure 4.17 where broken bound- 
ary lines represent Voronoi partitions when 
the distance is the drive time while solid 
boundary lines that when the distance is 
the spatial distance. 


Figure 4.17. Porcupine diagram super- 
imposed on a map. 


A similar approach to the above can be applied to help locate the Closed Circuit TeleVision 
cameras in a strategic location. Suppose that now the triangles in Figure 4.17 are important places, 
for examle business centres, schools, or universities; let us call these sensitive areas. Then security 
cameras should be placed at both points A and B as well as at the sensitive places. Here the solid 
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boundary lines meeting at the point A represent the Voronoi boundary for walk time instead of the 
actual spatial distance. Moreover, a camera should be placed at the point C which is the point where 
the shortest walkway linking two neighbouring sensitive areas meets with the Voronoi boundary of 
the walk time, in other words a point which lies midway along the shortest path that joins nextdoor 
sensitive areas. The reason for this is that these points are likely to be frequented by people going 
from one sensitive area to another, therefore they are prospectous for robberies. On top of that if 
a robber wants to flee from one sensitive area to another his sensible choice would pass through C 
while A with its connecting three different routes would be a clever choice to choose to confound a 
pursuer, or likewise B if he is driving. 

At a police headquarter and at the city council office I wanted to know whether security 
cameras have been placed at locations analogous to these points A, B, and C, but was told that the 
information is sensitive. The answer is understandable because at the time M. Bush was trying to 
get his hand on M. Bin Laden and M. Blair had just declared Britain at war and alert. However 
this was not good for my theory since we had no way of knowing whether this approach is new, 
or whether it could be used to improve the security in the city. Hopefully most of these strategic 
locations A, B and C have already been covered by security cameras even if the decision to do so 
had been arrived at by some other theories different from the one I suggest. Unlike in the case of 
sensitive places, the cost of the cameras to be placed at these strategic places may have to be the 
responsibility of the city council because their importance is not immediately obvious owing to their 
locations. In contrast, most schools, universities or businesses should be willing to help with the 
cost of such installation from which they would directly benefit. 

On one slide I explain how one can divide transformations into affine and non-affine ones. The 
affine transformations are governed by the equations x’ = az + by+p and y'’ = cx +dy+q. As 
a conclusion, mathematical models are beautiful because they give us understanding, and they are 
useful because they stem ideas and procedures necessary for solving problems. The presentation has 
gone from forest fire fighting to crime prevention without my knowing it, so the future works are to 
collect real data and practices, to analyse and compare with the theory and simulations, and then 
to implement and monitor the operations which I think unlikely to happen because it is impossible 
for me to become a member of the Greater Manchester police force. Dr. Jim Boran who was present 
gave good’s for project understanding, overview and structure, summary, examples, use of colour, 
connection with drawings, enthusiasm, pointing at screen, and familiarity with the material. He 
wrote that the explanations were clear, the content technical but with good personal view, but that 
I looked at the screen too often and at times was difficult to hear. Bigger labels are suggested. 

Naval forts are positioned with an idea similar to forest fire stations. Each fort forms a circle 
of effective defence around itself. These circles form a chain like beads on a necklace. One example 
(cf Petersburg, 2001) of this are the forts of Kronstadt near Kotlin Island in the Finnish Gulf, 20 
kilometres from Saint Petersburg, which have been used as a base for the Baltic Fleet and is now 
controlled by the Russian Navy. Here seventeen forts had been built between 1704 and 1896 which 
form a guarding barrier for the city of Saint Petersburg. Examples of these forts are the Totleben 
Fort, founded 1886-1913, and the Alexander I fort, founded 1836. 

In the study of traffic in § 8 we consider the percolation of roads, which are edges of the 
structure. We can also use the same structure for the study of fire protection and evacuation in 
cities, in cases of disaster, but now it is instead the percolation of cells that interests us. Roads form 
effective fire barriers as well as divide the conflagration into partitions. By carefully managing these, 
for example putting more efforts on the key or critical partitions, the flames can be contained more 
quickly and the damage minimised. These zones are also important for evacuation planning, for 
example that made by the Greater Manchester Police. In Tokyo where there is always a fear of big 
earthquake, informations are given out in plenty to tell people what to do in case such a disaster as 
the Great Kanto Earthquakewhich occurred on 1% September 1923 should struck. School grounds 
are usually assigned as evacuation centres. And as the distance one has to walk to get to one of 
these is crucial for one’s safety, the evacuation zones assigned by the city councils there naturally 
divide the city into domains that approximate the Voronoi tessellation. In 1998 when I was in Tokyo 
the fear reached its height, since it was generally believed that big earthquakes in Tokyo occur at a 
period of approximately seventy-five years. 


§ 4.12 Fractals 


In doing simulations on filtering membranes, one usually assume that the blocking particles are all 
of the same size. Or one may apply more intuition and experience, and say that particles should 
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have a normal distribution. Any mineral engineer will be able to tell that it is normally the case 
that the distribution of aggregates according to their sizes is normally normal. But one must keep 
in mind that the normal distribution here is by weight not number of particles. Therefore it is 
unlikely, except in a few special cases, that simulations which assume a normal number distribution 
of blocking particles will represent the real process. 

A pile of aggregates with a normal distribution is likely to have a fractal number distribution 
(cf Liebovitch and Scheurle, 2000). Such fractal distribution is hyperbolic towards the lower size 
ranges. But the overall distribution is more complicated, that is to say, a mixture of normal- and 
hyperbolic distribution. 

So far we have only considered percolation by competition between two phases. We have 
looked at networks in a cubic box and have assumed that the network is initially available in its 
entirety. But what if one has a competition of three phases or more? Or what happens when there 
are only parts of the network from the start? The picture of this can not be clearer than that of 
politicians arguing for votes, or the opinion polls of people across the social cross-section. Lustick 
and Miodownik (2000) consider this problem in the context of politics where the agreement clusters 
necessarily appear in various colours not merely black and white. 

Applying this to percolation and a whole new dimension opens up for investigation. One 
possible approach is to consider only an arbitrary part of the network within our cubic box. Clearly 
this is only meaningful if the part being considered is itself a percolating cluster. Carry this one 
step further and we can have percolation within percolation, when each cluster within the network is 
itself subject to another percolative process, threatening to destroy it for instance. Such a dynamic 
scenario is a yet untread water which should open up much ground for investigation. 

In the early days of the theory, literature on percolation only concerned itself with fluid wetting 
or blocking bonds or sites in the kind of situation that is said to be dual to the diffusion problem 
(Broadbent and Hammersley, 1957; Hammersley, 1957, 1961). The subject generally looked at self- 
avoiding walks and a measure of connectivity, the connective constant. From then on the appearance 
of the subject has somewhat changed, and it has become connected to Physics and fractals. The 
reason why there are very few works on percolation of a random network like the Voronoi tessellation 
is probably because there seems to be no needs for such study since the applications in Physics have 
told us that the regular lattices can represent the noncrystalline structure of the metal well. Moreover 
the only tools we have now, namely the Hamiltonian and power series method, can not be applied in 
its present form to random tessellations since there are no modulus relationships among the vertices. 
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§ 5. Porous media 

Patrik et al (1999) study the propagation in random Delaunay lattice. Particle on arriving at 
a site in the lattice deflects over the largest possible angle to either the right or the left, depending 
on the right- or the left scattering nature of the site. After the particle has passed through the site 
the latter goes into the reverse state. When a similar study is done on the triangular lattice, the 
entire trajectory quickly becomes confined to a particular strip which is bounded by two adjacent 
parallel lines of the lattice. They explain the propagation as being due to a blocking mechanism 
which prevents the particle from moving in a direction opposite to the propagation direction for 
more than a few steps. 

The flow of fluid in porous media follows the Darcy’s Law which states that the rate of flow 
through such a medium is proportional to the potential energy gradient within that fluid. The 
constant of proportionality is the hydraulic conductivity, which is a property of both the porous 
medium and the fluid moving through it, v = Q/A = —kdh/dl. 

The average velocity over the entire cross section is called the superficial velocity, vy = Q/A = 
tn/pA. The velocity that is based on the actual open space within the porous media is called the 
interstitial velocity, vn = Q/(e€A) = m/epA, where ¢ is the porosity, e = V,/Vs = (Vi — Vs)/Vi, Va, 
V, and V; are respectively the void, solid and total volume. 

The definition of porous media is different from that of porous materials. A porous medium is 
a medium through which other substance may pass, whereas a porous material merely means some 
certain kind of material the internal structure of which is filled with pores. Most of porous materials 
can be used as porous media. All of them are useful because of their internal structure, and though 
we do not need to know this to be able to use them, we do have to understand it if we want to use 
them efficiently, or if we want to improve upon some of their particular properties. 

The outer bark of Quercus suber L., for example, has a property that is ideal for its commercial 
use as cork. Its material property is such that it neither shrink nor expand when stretched or 
compressed. This makes it an excellent insulator of both heat and sound as well as good for damping 
vibration (cf Ashby, 1990). It has high frictional coefficient, and is both impervious to liquids and 
chemically stable. These properties make it ideal as wine corks. By understanding its internal 
structure and mechanism that makes it stay in the same shape when experiencing external forces, 
scientists have succeeded in making a synthetic material that shrink, instead of expand radially when 
pressed. When such material is used as cork, it can be easily pushed inside the neck of a bottle. 
When we release the force it will expand to its full size and fit snuggly where we placed it. 

Pores introduced into a fabric can also increase its commercial value by giving it a lighter 
weight, higher heat retention rate, which makes it ideal for both ski and snow-board wearer and 
casual wear (JETRO, 2001). 


§ 5.1 Zeolites 


Zeolite is the name given to a group of minerals with a porous structure. It is used as molecular sieves 
as well as in many chemical engineering processes. It may be an important ingredient responsible 
in lubricating the fault lines in earth crusts, thus making earthquakes less severe. Evan et al (1995) 
found foliated cataclasite and ultracataclasite in the San Andreas faults. These are composed of clay 
and zeolite. In the case of the ultracataclasite, there are fragments of 20-100 ym diameter feldspar 
and quartz embedded in a matrix of clay and zeolite which has grains of sizes smaller than 10 wm. 

The chemical composition of zeolites is that of sand, ¢.e. aluminosilicate. Their general formula 
is (Dana and Dana, 1997) 

(Na, Ca, Ba)(1-2) (Al, Si)5O10 5 nH2O, 
or, as given by Gottardi and Galli (1985), 
(Li, Na, K) (Mg, Ca, Sr, Ba)q[Al(a42a)Sin—(a4+24)O2n] -mH20. 

Normally m < n. They differ from sand in that they have large internal cavities. These cavities 
increase their internal surface area as a result of which they can be used as catalyst or molecular 
sieves with both shape- and size-selectivity. One example is the use of zeolite as sieve to separate 
iso octane, which has high antiknock property, and iso pentane from octane and pentane. Zeolites 
can not withstand high temperature because they rapidly lose water. Even though this process is 
reversible, their molecules collapse at temperature higher than 600°C. Around 41 types of zeolite 
occur in nature. Linde Division of the Union Carbide Corporation produced the first synthetic 
zeolites in 1950s. 


152 


Ph.D. Thesis, UMIST. K N Tiyapan. Chapter 5: Porous media 


The name ‘zeolite’ was coined by A. F. Cronstedt in 1756 from the Greek zein, to boil, and 
lithos, stone. The structural units of zeolite comprises of the primary building units of silicates in the 
form of XO, tetrahedra, where X is mainly Si. Other basic units are the chain of fibrous zeolites, the 
singly connected 4-ring chain, the doubly connected 4-ring chain, the 6-ring either single or double, 
the hexagonal sheet with handles and the heulandite unit. 


§ 5.2 Crystalisation 


In the diffusion theory of crystal growth, the overall rate of crystallisation is determined by 
the rates of two processes occurring one after the other, that is diffusion of solute from the bulk 
solution to the interface between solution and the crystal, followed by integration of solute atoms 
into the crystal lattice. The rate of diffusion of the former is Ry, = ki,(c — cj), while the rate of 
surface reaction or integration is R, = k,(c; — c')", where c; is the interfacial concentration. The 
value of the diffusion mass transfer coefficients k!; when estimated from the growth is considerably 
different from the same value that is obtained from the dissolution experiment. But if we assumed 
that they are equal, then R, = k,(Ac — R,/ka)”, which, when n = 1 gives R, = K Ac” where 
1/K =1/kg +1/k,, and when n = 2, Ry = ka [(1 + ka/(2h,-Ac)) — [(1 — ka/(2kpAc)?) — 1]1/?] Ac., 
but there is no general solution for all n’s (cf Garside and Mullin, 1968). 

Particles in a crystalliser are kept in suspension by a stirrer or a pump. Lim et al (1999) 
study the effect of the attrition when crystals encounter with high speed impellers inside these 
apparatuses. When the impact energy exceeds the crystal strength, crystals fracture which gives 
rise to a particle size distribution as a result. Impact can occur at faces, edges or corners of crystals, 
but in the attrition model only the contact of a crystal corner with another flat, much harder object, 
for instance the steel impeller, was considered. Furthermore, crystals are assumed to have the shape 
of a cube or other polyhedra, but in the model the crystal faces forming the contacting corner are 
replaced by a cone having an included angle of 120°. In the vicinity of impellers the flow is turbulent 
and crystals travel in a random manner. Repeated attritions reduce the crystal size, assuming no 
competing effect of crystallisation. They assume a normal volume shape factor distribution and 
proceed to simulate according to a procedure recaptured here as Algorithm 5.1. The input to the 
algorithm is the impact energy W,; the outputs are wy(J) and c(J). The hardness of a solid is 
its resistance to local plastic deformation. The contact pressure of a plastically deformed cone is 
assumed to be the same as the Vickers hardness H,. An isotropic material has two independent 
elastic constants, namely the shear modulus p and the Poisson’s ratio v. All the other constants can 
be determined from these two quantities. Here a is the volume shape factor, I the fracture resistance, 
K,, the efficiency of stress field created by the impact between the crystal and the impeller, « stress 
field parameter, & = 5, ty quasi-isotropic shear modulus, N the total number of fragments, W, 
impact energy, H dynamic hardness, Lmin and Lmax respectively the minimum and maximum size 
of a fragment, a characteristic size of the plastic zone, r the distance from the peak of the cone, rmax 
the maximum distance from the peak of the cone to the newly created surface and J the particle 
size classes. The values of Hy, wy» and W. used were that of magnesium sulphate heptahydrate and 
potash alum; r; and rg are random numbers. 


Algorithm 5.1 Fracture by attrition, Lim et al (1999). 


for 1 = 1 to 50 do 

ater; 

T/K, © «We? H®/3/(5.2u) 

N «7x 10-*W,H°K3/(ap*T?); 

for j = 1to N do 
Imin + 32uT'/(3K,H?); 
Emax “ Tmax /23 
calculate rmax} 
r + exp[(13 log a — log ra) /13]; 
L © 3plrt /(WiP HH? Kr); 
find J corresponding to L; 
wewt+L; 
wi) — w(J) + LD; 
c(J) H cS) +1; 


endfor 
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endfor 


wy(J) © w(J)/w3 : 


§ 5.3 Fluid flow within networks 


Flow of viscous fluids through networks of geometrical objects is an important topic in various 
disciplines of engineering. Happel (1959) studies two cases of the flow of viscous fluid relative to 
arrays of cylinders, one parallel while the other perpendicular to the cylinders. This is found in 
practice as the flow through a bundle of heat exchanger tubes. 
The flow pattern within a void of a porous media follows the equations of flows with vorticity, 
in 2-d, 
Oy  1d% Ea er, 
Or? or Or or? O02” 
and in the axisymmetric case, 
QO Ow “ O OD | 
ror sinOOr  r300sin@d0 
Or equivalently to both equations, 
és 1 (2 = a) 


r Or 06 


when vg = 0%/Or and v, = —Oy/(rd0) in 2-d, and vg = OW/(rsin@dr) and vp = —Ow/(r? sin 008) 
in the axisymmetric case. The vorticity, ¢, can be used as boundary conditions, for instance ¢|-=9 = 
—0v,/a00 or C|r—a = Ov, /a00. Or it can be used within the hole, for example ¢ = k sin @ which gives 
the solution w = (eyr7! + egr + car?) sin @, or = = (egr + car?) sin @ sin @ if c, = 0. This solution is 
a part of the general polynomial for the stream function, 


4 (cf Rowe, 1965.) 


pa=(r-cor *+ar +e +egr + cur? +--+) sind. 


§ 5.4 Material science 


A chiral is a group of points, or geometrical figure, whose mirror image can not be brought to 
coincide with itself. Achiral is the antonym of chiral. A chiral object is an object which fails to be 
achiral. Chirallity is a purely geometrical property since since all the operations involved, namely 
the plane symmetry of the reflection, the rotation, and the translation, are all isometries. Chiral is 
used in the studies of molecules and knots. 

The Poisson’s ratio, v, is the ratio between the transverse contraction strain and the longitudi- 
nal extension strain, that is vy = —e,/<,;. The theoretical value of the Poisson’s ratio can range from 
-1 to 0.5, but for normal materials it is generally positive. At v = 0.5 the bulk modulus is much 
greater than the shear modulus and the material is incompressible, while at v = —1 the opposite is 
true and the material is very tough and highly compressible. The bulk modulus B, and the shear 
modulus G are related to each other by the equation B = 2G(1 + v)/(1 — 2v). For rubber this 
value is 0.5, for aluminium 0.33, while for cork it is approximately zero. Materials with a negative 
Poisson’s ratio have been found whose structure is re-entrant (Lakes, 1987). Applications of such 
materials include robust shock absorbing material, fasteners, and stoppers of the wine bottles. 

Chemical reactions can be thought of as phase changes in percolation. For example the ex- 
traction curves for Sb shown in Tiyapan (2003, KNTs(ii), § E.1). These curves represent an sshape 
starting from one phase, represented by 0% Sb, to another at the maximum per cent extraction 
where it saturates. All the per cent extraction graphs shown in Tiyapan (2003, KNTs(iii), § E.1) show 
them with high slope at t = 0. This is only for the convenience of drawing, since it is difficult to 
draw curves with an s shape smoothly on the computer. Also, the data obtained from the experi- 
ments do not extend to the time immediately following ¢ = 0. Extracting the solution using pipette 
takes some time to do, and therefore it has not been possible to prove experimentally whether the 
extraction yields start off with zero slope at ¢ = 0 or not, even though one might conjecture that 
this is likely to be the case if one considers the extraction yield in the leach solution as a developing 
phase in the continuum of the solution. 

Particle sizes come in a variety of definitions. Svarovsky (1977) divides them into three groups, 
namely definitions by equivalent sphere diameters, by equivalent circle diameters and by statistical 
diameters. These are listed together in Table 5.1. 
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diameter 
volume diameter 
surface diameter 
surface volume diameter 
drag diameter 
free-falling diameter 
Stoke’s diameter 
sieve diameter 
projected area diameter 
projected area diameter 
perimeter diameter 
Feret’s diameter 
Martin’s diameter 
shear diameter 
maximum chord diameter 
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criterion 
equivalent volume of sphere 
equivalent surface of sphere 
surface to volume of sphere 
resistance to motion of sphere in the same fluid at the same velocity 
free-falling speed of sphere, same fluid and particle density 
free-falling speed of sphere if Stoke’s Law is used (Re < 0.2) 
diameter of sphere passing through the same square aperture 
projected area of a circle, the particle resting in a stable position 
projected area of a circle if the particle is randomly oriented 
perimeter of the outline of a circle 
diatance between two tangents on opposite sides of the particle 
length of the line which bisects the image of the particle 
particle width obtained using an image shearing eyepiece 
maximum length of a line limited by the contour of the particle 


Table 5.1 Particle size definitions. 

Particle size distribution comes in four types, distribution by number f,,(x), by length f;(x), by 
surface f,(a) and by mass or volume f,,(x). Conversions among them are done by f; (x) = kitfn(a), 
fs(a) = koa? fn(x) or f(x) = kg? fy(x). Conversion is only possible if we know the shape factor’s 
dependence on particle size, because k; often contain a shape factor. The distribution frequency is 
by definition {5° f(a)da = 1. 


(m,df(a)/dx,Q) +224 


overflow 
—> 


—> 


separator (my, df¢/dx, O) 


| underflow 
(mc, df. /dz,U) 


Figure 5.1 Schematic diagram of a separator. 


Small particles tend to flocculate with one another. This makes it impossible to do experiments 
with very fine powder. Therefore any experiment which claims too fine a size as a control parameter, 
for instance in the micro metre range, must bear in mind that particles of such small size ranges 
tend to cluster into hard agglomerates bound by strong chemical bonds, for example from a previous 
chemical or thermal treatment, or into soft agglomerates by van der Waals attraction or by capillary 
forces. In ceramic making, this results in nonuniform packing during forming which leads to large 
voids or flaws after thermal processing. As one expert put it, you can not have powders in the size 
range of microns because they will disappear instantly into the air. One way to deal with them is 
to mix them into slurry. 

Powder is also subject to impurities from the earliest stage of its life, that is during the 
comminution, which can be done by various methods, for example mortar and pestle, ball milling, 
jaw crushers or crushing rollers. The impurity comes from abrasion of the grinding media, which 
can be effectively avoided in the case of the jet mill, where particles are driven by air streams in 
opposing directions to collide with one another among themselves. Powder can be characterised by 
its size distribution, surface area, shape, composition and crystal structure. 

Particle size distribution can be determined by means of Stoke’s law, screening and microscopy. 
Stoke’s law states that the steady state velocity of a spherical particle travelling through a fluid is 
Ur = 2(Pp — pp)ar?/9n, where a can be due to gravity, in which case a = 9.8 m?/s, or alternately 
due to centrifugal force, a = w?x. The time required for a spherical particle of of radius r to travel a 
depth x is t, = z/v,. No particles of radius larger than r will be found at depths less than x at time 
t,. The fraction extracted mass at short time intervals in an experiment using a graduated cylinder 
is constant, 7R? fins fodx =m. 

The equivalent spherical radius is the radius of a solid sphere that has the same steady state 
velocity in the slurry as the particle. The volume distribution function is f,(r) = dV/(Vrdr). 
The number of particles in the powder having radius r is N(r) = 3Vrf,(r)/4ar?, the mode of the 
distribution is the peak and the mean radius is [5° rN(r)dr/ f)° N(r)dr. 

When finding particle size distribution using screening methods, sieves are arranged in a stack 
with their apertures decreasing from top to bottom. Mass fractions retained on each of the screens 
is f; = m,;/mr, which can then be converted to N(r). 

Microscopy methods use a variety of microscopes, for instance optical, SEM or TEM, depending 
on the magnification required. Measuring the dimensions and determining the shapes of the particles 
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can be done manually with a ruler, or by using an image analyser. 


Powder surface area is determined by gas adsorption methods. The specific surface area, i.e. 
surface area per unit mass, can be found by s = NoV,,/V., where N is Avogadro’s number, oo the 
effective cross sectional area of the adsorbate molecule, V,, the volume per gram of a monolayer at 
STP and Vo the STP number, 2.24 x 104 cm?/mole. 


Monolayer adsorption is used when molecules strongly interact with a surface to in a uniform 
monolayer. The Langmuir equation is 0 = kP/(1+kP), where 0 = V/V,, is the fraction of sur- 
face covered and K a constant. Then P/V = 1/(KVm) + P/Vm and the Langmuir isoterms are 
plots between V and P/F, Po being the saturated vapour pressure. For multilayer adsorption, 
“/ [V1 —2x)| =1/(eVn) + #(e — 1)/cVm, where & = P/Pp. 

Permeametry measures the flow of gas under a pressure head through a packed bed of particles. 
The Darcy’s law of laminar flow through porous medium is u = (1/A)(dV/dt) = KAP/€. The 
Poiseuille equation for a fluid flowing through a tube of radius r is @ = r*AP/(8n£), where 7 is the 
viscosity. It is also written as h/L = 32nv/(gd”), where h/L is the headloss per unit length, 7 the 
kinematic viscosity (m?s~+), 7 = /Pw, Pw mass density of the water (kg-m~3), y absolute viscosity 
of water (kg(ms)~!) and d the diameter of the particles. The hydraulic radius is the ratio between 
the wetted area and wetted perimeter, r = (md*/4)/(ad) = d/4. Therefore, h/L = 2nv/(gr?). 
The number of particles is n = Vy/V;, and the total wetted surface area of solids is Ap = nAj. 
Therefore Ap = (1 — €)(ad”)/(md?/6) = 6(1 — €)/d. Then, replacing the factor 6 with the shape 
factor s, r = ed/[s(1 — €)]. Because v,A, = vy Ay/R, where p and f means pore and filter surface 
respectively, R is a constant which accounts for friction loss, and also J = 1/R, therefore we have 
h/L = 2nvp/(gr?) = In(l —€)? vps? /(gerd?). 

Kozeny used the ratio V/A = d/4 for a cylinder and the porosity p = V,/(V, + Vs) and 
gavethe equivalent cylindrical diameter for a packed bed as d. = 4pV,/[A(1 — p)]. Then the average 
velocity through the channels is up = dvp/dt = [(V2AP)/(2A?nl)] (p/(1 — p))?, where | is the 
channel length. The Carman-Kozeny equation accounts for the tortuosity of the channels in the 
bed, s? = 2.5p3AP/(2pnlv(1 — p)”), where V; = m/p, the specific area s = A/m and p the solid 
density of the powder. The determination of the phases and composition of a powder is done by 
x-ray powder diffraction, thermogravimetric analysis and differential thermal analysis. 


Wet methods are preferred when preparing material before firing in order to control agglomer- 
ation during forming, for otherwise when the particles are agglomerated the narrowest particle size 
distribution becomes of little or no help. In wet forming, solid particles are suspended within a liquid 
the chemistry of which can be adjusted to make them mutually repulsive or mutually attractive, 
for respectively deflocculation and flocculation. Stabilising makes particles mutually repulsive either 
by adsorbing polymer chains on to the particles in steric stabilisation or by putting charged ions 
or polar molecules on the particle surface in electrostatic stabilisation. In steric stabilisation, one 
end of the long polymer chains, which is hydro- or lyophobic, is adsorbed on the particle while the 
remaining end, which is hydro- or lyophilic, extends in to the liquid. 


Electrostatic stabilisation results in diffusive double layer built around particles. Zeta potential 
measures the repulsive potential between particles travelling in a medium. It defines the electrical 
potential in the double layer at the surface of shear slippage when the particle is forced to travel 
through the fluid when an electric field is applied. The composition and thickness of the double 
layer can be changed by adjusting the pH of the suspension. The point where the pH gives a zero 
zeta potential is called the isoelectric point. It is a point where spontaneous agglomeration occurs. 
The electrophoretic mobility vy. = u(E)/E = Cepeo/(fnn), where u is a steady state velocity, fr a 
constant determined by the dispersing medium and the particle size and 1 < f, < 3/2. 
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§ 5.5 Forces between particles 


Essential to the simulation of stochastic models of particles is the consideration of forces acting 
on each particle (cf Schumacher, 1996). Since in theory the forces acting between two particles 
extend their effects to infinity, one has to make approximations. The degree of justification to these 
approximations depends on how the force in question vary with distance. 


10” 


In the case of the Coulomb force where the 
jee ER ALTE SOS eee oe “=... | potential energy is u(r) = qi qo/(4me0err) and 
= the force f(r) = —du(r)/dr is 


5 f(r) = q1¢2/(4r e067”), 


the justification is high as is seen in Figure 
5.2 where the force-distance is a decreasing 
straight line on the log-log graph, qi and qe 
are protons, and r the distance in metre. 


Figure 5.2 Coulomb force vs distance in log- 
log scale.. 
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The effects of the Coulomb force are prominent when sizes of the particles are small. In 
quantum mechanics where the scale is atomic, for instance, the problem becomes one in which there 
are a large number of particles interacting with one another, i.e. an n-particle problem. In theory 
the Hamiltonian operator can be applied and then the problem solved numerically. But in general 
this is not possible due to the too many particles involved in the calculation. Theoretical solution 
is possible by the various methods of approximation, for example the Hartree method (cf Brown, 
1972) where the best wave function is found in terms of the one electron function, i.e. orbitals, ©, 
or a Hatree-Fock approximation which reduces the number of equations to n/2. 


§ 5.6 Arbitrarily shaped particles 


An arbitrarily shaped particle in a divided space will necessarily have a harder time travelling 
around compared with a sphere even if, or rather especially when the space is a mathematical one. 
Small particles are roughly spherical, or so the lore of science says, but when one’s technology has 
led one down to the realm of dimensions in the order of those small particles, for instance the nano 
technology, then shape does matter a great deal and the universal spherical assumption can no longer 
suffice. 


This is only one of the reasons which justify the investigation into the cases where the ratio 
of surface area to volume is not a minimum. When physics and motion are involved, this kind of 
study can incorporate both the continuous trajectory calculation as well as the discrete raster grids 
of the digital visualisation technology. The best way to hold such data is as matrix maps or masks. 


The best and quickest way to define the boundary of an irregularly shaped object while allowing 
ease of investigation of its interaction with the surrounding is as raster outlines. To do this, I first 
investigated such outlines computed from the equation of straight line segments which make up the 
boundary of the object. Hereafter object boundaries are defined without loss of generality as straight 
line segments that link all the vertices together with no gaps. The first investigation is unsuccessful 
as it violates this definition by leaving many gaps, as can be seen in Figure 5.3 (a). The second 
algorithm calculates the grids from the coordinate axis along the direction of which the slope of the 
edge is minimum. But this too still leaves gaps along the boundaries. 
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Oe “2 F Hf i 
(a) (b) 
The object shape is defined by a set of landmark j i 

vertices or corners. This must be listed in intended fe - ra 

order, otherwise it would be impossible to deduce 
the order intended unless one assumed that the ob- 
ject is convex. Figure 5.3 is an example of irregular — 
object shapes, the name Kit in Anglo-Saxon’s runic 
script (cf Freeborn, 1992) created by three different f 4 
algorithms. j 


= ib, H 
(c) i (| : 
Figure 5.3 Kit the existential object. 


There are three different programs in § A.24. Only the third one, tioa3.m gives a satisfactory 
result which leave no holes in the surface of the objects. In the second program, tioa2.m, v is a 
data set which contains the number of vertices, the vertices, slope, map matrix and the object’s 
dimension. In tiao3.m, v is somewhat different. It contains the number of vertices, the vertices in 
real dimension, the grid dimension of the box containing the object, the vertices in grid dimension 
and the matrix map or image of the object. In order to exploit possible parallelisation, the last 
program ignores the warnings when there is a division by zero as this means that all the grid 
positions of the whole line can be operated upon together as a vector quantity. 

The third program operates as Algorithm 5.2 does. It utilises the same idea as that used by 
a child or an artist alike when they draw or paint. Drawing and painting are both one-dimensional 
process which seeks to produce two-dimensional results. Each pencil- or brush stroke travels along a 
path the direction of which has one dimension. Unless he uses a very thick brush, there will always 
be a possibility of gaps forming between stroke lines, which can be annoying because however small 
they may be one needs to use a great deal of paint in order just to cover them up. This problem can 
be overcome by painting along two directions, for instance perpendicular to each other especially 
when the person concerned is a child. At a first glance, or to a novice, this may seem a wasteful 
practice, especially with dear types of ink. But a little practice and experience will show that this 
proves in most cases to be more economic than painting along only one direction. And since it can 
be readily seen that painting in two directions perhaps already uses twice the amount of ink required 
to cover the paper, it follows that doing so in one direction only would use a great deal more ink 
than this. 


Algorithm 5.2 Particles or objects with arbitrary shape. 


(i,7) < find grid coordinates for (x,y); 
find all vectors linking vertices; 
{pi} discretise these vectors into sets of points; 
for all intervals between consecutive p;’s do 

draw along y-direction to fill gaps; 

draw along x-direction to fill gaps; 
endfor 

o 


The problem of texture quantification is in the difficulty in expressing such quantity by a single 
parameter. Properties such as coarseness, smoothness, heterogeneity and regularity are finger prints 
of particles, and are the result of physical and chemical processes. By using fractal geometry, it is 
possible to describe quantities like texture. A fractal is a set whose metric properties can only be 
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consistently illustrated with a dimension D > T, where T is the standard topological dimension. 
We can express this dimension as D = T + (1 — H), where H is the codimension. Curves can have 
their roughness described. By giving a fractal number between one and two the space filling ability 
of the curves is established (cf Van Put et al, 1994). The first method defines the fractal dimension 
of a function in terms of the Fourier power spectrum P(w) = w7@#t)); H and D are determined 
by doing a linear regression on the log-log plot of the observed power spectrum as a function of 
frequency. Another method defines fractal dimension in terms of how the variance of interpixel 
differences changes with distance; D is estimated from a log-log plot, a variogram, of the variance 
of increments versus increments o?(x) & x?/. 

Techniques used for measuring sizes of particles include sieving, microscopic analysis, electronic 
particle counters, laser diffraction analysis, permeability-, sedimentation and elutriation methods. 
Usually arbitrarily shaped particles are characterised by a variety of methods. For a single particle 
this usually means transforming its property to the corresponding value of a sphere. A particle may 
be represented by a sphere which has the same volume, surface area, surface area per unit volume, 
or the area projected perpendicular to the flow direction. It may be compared with a sphere which 
has the same settling velocity in the same fluid, or a sphere which will just pass through a square 
aperture of the same size. Or its area projected on to the ground, when it is resting in the position of 
maximum stability, may be compared with that of a sphere in the same situation. To be consistent 
in one’s choice is more important than which choice one chooses from. The sphericity of a particle 
is defined as © = Aj,/Ap, where Aj, is the surface area of sphere which has the same volume as the 
particle and A, surface area of the particle itself (cf Coulson et al, 1991). 


§ 5.7 Non Poisson number distributions of particles 


Let us look at an aggregate of spherical particles whose average diameter is 200 + 50 ym, 
density p = 3.7, and the total weight one gramme. Assuming the gravitational constant to be 9.81 
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The overall weight is W = n-w, where w = pVg 
is the weight of each particle. The particles be- 
ing spherical their individual volume is thus V = 
(4/3)ar? and their distributions by weight and by 
number are shown in Figure 5.4. Figure 5.4 shows 
(a) distributions by weight, (b) distribution by num- 
ot | ber in normal scale and (c) with log scale in the 
y-axis. 


a 


number of particles 


Po 
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Seer ao Figure 5.4 Distributions by weight and by number. 


In particle size classification the separation size, D509, is the size where there is an equal chance of 
particles being reported to the fine or coarse fraction. The grade efficiency curve is normally plotted 
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between the weight fraction to coarse product and the normalised particle size D/Ds0. The choice of 
classification equipments depends on many factors, for instance the size of the particles in question 
and their electrical or magnetical property. For fine particles, Treasure (1965) discusses several 
types of classifier, namely the solid-bowl centrifuge, the Hosokawa micro-separator and the Head, 
Wrightson air classifier all three of which have Do « \/¢/w, the hydrocyclone where Dso « 1/ g®® and 
q x p®, and the Alpine Mikroplex classifier where Dsq « u/(wVV) where u and w are respectively 
the radial and tangential velocities and V is the volume air flow rate, and Dso « 1/w in practice. 


The distribution of bubble size in gas fluidised beds had been represented by various kinds of 
distribution, for example the log normal, Gamma and x? distributions. Rowe and Yacono (1975) 
preferred the gamma distribution of the volume, that is 1(V|m,n) = [1/(n™I'(m))] V™—1 exp(V/n), 
in its normalised form I'(v) = [m™ exp(—mV)V™—!] /T(m). 


Sieve testing is affected to blinding materials, that is parts of materials under test which lodge 
themselves in the apertures of the sieve. To account for this Rose and English (1973) give a general 
rule of thumb which says that a particle for which @ < tan~!y will blind the sieve. Here @ is the 
angle that the line of reactive force exerted by the material on the particle, which is normal to the 
contact surface, makes with the horizontal plane. Then they give the empirical value of 7 to be 
approximately 20°, which results in 6 < 0.3358 radians as being the criterion for blinding to occur. 
For a spherical particle this gives the ratio x/r = cos @ > 0.9441, where r is the particle’s radius and 
x half the minimum distance between two normal reactions which act on the particle. Therefore the 
radius of particles which blind the sieve must be such that r < 1.0592x, which corresponds to the 
aperture of 1.0592 times greater than the aperture if we assume that the aperture is 27 through out, 
or equivalently that the opening has its corresponding faces vertical. 


It would be of a theoretical interest to prove this equation 9 < tan~!y given by them. Also, 
the angle of friction, 7, is not a constant but varies for different pairs of materials in contact with 
each other, as well as changes the value when we change the medium surrounding them to a different 
type, or even when the moisture content changes for that matter. All of this goes to say that the 
value 1.ldmin extensively used in chemical engineering literature is probably nothing more than an 
engineering rule of thumb based on two conjectures, one in the limiting value of 6, the other in the 
value of y. But it gives a simple and convenient criterion for use when doing computer simulations, 
which are nothing but means to visualise the conceptualised physics anyway. Other values of @max 
should give the same qualitative result. We may, for instance, choose the coefficient of our dmax to 
be 1.02, corresponding to 6 > 1/cos10°, instead of 1.06 which corresponds to the 9 > 1/ cos 20° 
used above. 


Since no explanation of the formula regarding @ given by Rose and English has been given 
elsewhere, I arrive at my own derivation as follows. First, let us recall from our sixth form physics 
how friction can be described by either the coefficient of friction or the angle of friction, both of 
which are merely different sides of the same coin. The angle of friction is also known as the angle of 
frictional resistance, the internal angle of friction and the angle of shearing resistance. For clean sand 
it is approximately the angle of repose. It reduces with moisture content and is zero for saturated 
clay. For solid it is approximately the angle of inclination of the surface of one material at the point 
when a block of another material, placed on top of the former, starts to slide. This is precisely 
the method by which the coefficient of friction is determined. There are two different coefficients of 
friction for each pair of material, namely the static- and the kinetic coefficients of friction, the former 
having a slightly higher value than the latter, and both are defined as being the ratio between the 
limiting friction and the normal reaction, p = f;,/Fn. This gives rise to the formula for the friction 
being f, = uF’, which is equal to mg sin 6 at the point when the object of the second material slips, 
0 being the inclination of the plane made by the first material. Therefore we have = tan 6, where 
ps is the coefficient-, and @ the angle of friction. 


Coming back to our formula in question, 6 < tan7!y. Figure 5.5 shows a spherical object 
sitting on top of an aperture in a sieve. From Figure 5.5, a+ 6 = 1/2. Without loss of generality, 
assume 7 < 7/4. Let up = tana, be the coefficient of friction at some critical value a,. Furthermore 
assume that a, is small, so that we can approximate ys by the angle of friction y. Then we have 
a, = tan-+y = tan! y. It follows that the particle will block the pore whenever a < a,, whereas 
it will pass through the latter when a = 7/4. 
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Because we know from experiments that some of 
the particles does lodge against, and blind the pore, 
BeOS SS there must exist some a, < 8 < 7/4 such that this 
’ % occurs. For the reason that the number of parti- 
i \ cles that blind is empirically small we know that 
i \ B is small. Since we already have the relationship 
among a, @ and 7/4, namely a+ = 7/4, and since 
‘ ' Q- is also small from our assumption, then in or- 
, a) ef der not to unnecessarily introduce another param- 
aera Sag DS alate eter into what necessarily already contains some 
i error due to all the approximations so far made, 
we let 6 = a-. Now we may say that the parti- 
cles will blind if and only if 90° -a, = 90° -6B< 
B, a < 90°. For other remaining values of a, viz. 
\ both a, < a < 90° and a < a, in other words 
‘ 0<a< 90° —a,, then, particles must block. But 
\ 9 = 90° — a, therefore the particles pass through 
Oh. ._ ‘ when 6 = 0, blind when 0 < @ < a¢ and block 
Ss ‘ when a, < 6 < 90°. So the blinding particles 
—s.. 0(a) have @ < a, = tan7!¥y and the approximation is 
a.(0) ee x explained though not q.e.d.’ed. 


Figure 5.5 Spherical object on a sieve. 


Let A, be the available area free from blinding material at some instant, A the physical area of 
the sieve cloth, G the area of cloth blinded by a unit mass of blinding material, and wo and w the mass 
of blinding material on the sieve respectively at t = 0 and at some t > 0. Then A, = A—G(wo —w) 
and dw/dt = —kiwAg, where ky is the diffusion coefficient, and then dw/[w(b + w)] = —adt, where 
a = —k,G and b = (A/G) — wo. Assume b 4 0. At t= 0, w = wo. It follows that 
= bwo exp(—abt) 
~ b+ wo(1 — exp(—abt)) 

Furthermore, let W be the total weight of material on the sieve, W. the residue on the sieve 
for a theoretically infinite time of sieving, and K a constant. Then, dW/dt = —K(W —W,.)Aa, and 


it follows that dW/(W — W..) = —KG(b+ w)dt = —KGbdt — KGbexp(—abt)dt/(q — exp(—abt)), 
where Q = (b/wo) +1. Again assume b 4 0, and W = Wo at t = 0. Then, 


In [(Wo — Woo )(W — Wo) =K {Gbe+ (E) py Coen ah (22). 


and then W = W.. + (Wo — Woo) d, where X = b/wo = (A/Gwo) — 1 and 


fs exp(arwot)(A+1) 1 =HGiP 
> Xr a : 


When the total amount of blinding material is sufficient to bind the sieve completely, A/Gwo < 1 
and so —1 < X < 0, when it is exactly enough to completely glind the sieve, A = 0, otherwise 
A/Gwo > 1 and 0 < A < o, Gwo being the area of sieve cloth which will be blinded by all the 
blinding material present. At t + oo, 5 = 0 when 0\ < oo and 6 = (—A)*&/*, a positive fraction, 
when -—1<A<0. 


w (21), 


(23), 
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§ 6. Filtering membranes 

Tanemura et al (1983) give a good algorithm which is both powerful and simple for constructing 
the three-dimensional Voronoi tessellation. The circumcentre of a DT being a vertex of the VT and 
the circumsphere of a DT being empty, this second condition of which is known as the contiguity 
condition, they arrive at an algorithm which works for both nondegenerate and degenerate cases. 
Algorithm 6.1 summarises their algorithm. Associated with the atom i, let V; be the polyhedron, 
vu, vertex atoms of V;, and S; set of atoms surrounding 7, $; D v;. Here H;(a6|7) is the half-space 
determined by {i,%,,ig} and does not contain 7,, S;(aB|y) C Si C Hi. 


Algorithm 6.1 Voronoi construction in three dimensions, Tanemura et al (1983) 


find 7,, nearest, and thus contiguous, to 2; 
Cee a3 
find j = ig s.t. {i,i1, 72} has the minimum circumradius among all {i, 71,7}, 7 € Si, J # 41, 123 
Cy 2g UC 
find j = i3 s.t. {i,%1,%2,i3} has the minimum circumradius among all {i,i1,i2, 7}, 7 € Si, 
GF ti, ta} 
Vi = {i, tr, ta, t3}5 
Cy — ig UCi3 
clear {m.} and {lag}; 
al; 
for all ig € C; do 
ifm, = 1 then 
atat+l; 
else 
find {i,i9,ig,t,} € T;, which has {¢,i,} in common and where [gg = 1; 
SilaBly) — SiN Hi(aly); 
describe a circumsphere to each quartette {#,%.,%8,j}, where j € Si(af|y) and m; #1; 
find jmin for which the centre of the circumsphere of {i,49,%3, jmin} has the minimum 
signed z-coordinate value among all circumspheres obtained, z-axis being normal to 
{i,iq,i8} and lies towards the side of H;(af|y); 
T; + {i, ig, if, is} UTs 
if is ¢ C; do 
Cy — 16 U CS lag + +3 las +43 Ugg + +3 lag +45 Usa + 4; lag + 4+; 
endif 
if all {J,} are equal to 2, mz, + 1; 
if all {Ig} and/or {ls} are equal to 2, mg + 1 and/or ms + 1, 
endif 
endfor 
find geometrical quantities of IL,; o 


Algorithm 6.1 is justified by proofs of theorems which confirm that atoms nearest to each 
other are contiguous to each other, a triangle with the minimum circumradius is a face of a DT, 
which means that all its three atoms are also mutually contiguous, likewise a tetrahedron with the 
minimum circumradius is a DT, all its four atoms are mutually contiguous, two DT’s sharing a plane 
have their fourth vertices at a signed minimum distance on either side of the plane they share. Some 
interesting geometrical lemmas help towards the proof of these theorems. 

Jackson (1994) studies porous media and represent them using Voronoi models. His interests 
are in the synthetic membranes, which can be symmetric or asymmetric, homogeneous or heteroge- 
neous structure, neutral or charged, and passive or active transport. The resistance to mass transfer 
is mainly in the top layer if it exists. This is a dense layer 0.1-0.5 wm thick which lies on top 
of a porous sublayer 50-150 ym thick. Regular tessellations in two dimensions come in 11 distinct 
tilings called Archimedean tilings. Tessellations by random lines in 2-d have the equation of the lines 
xcos6+ysin@—h = 0, where 0 < 6 < m and —o0 < h < ~w, which is analogous to tessellations of ran- 
dom planes in 3-d with the equation of a plane x sin cos¢+y sin é sin é+zcos@ =0, where0< ¢< 
a. The definition he gives of the Voronoi polygon is (es HA («;,2;), where the half-plane H (xj, x;) is 
{x € E?|d(x,x;) < d(z,2;);i #7}. The bisector is B(a;,x;) = {x € E”|d(x,x;) = d(x,2;)}. There 
are n. + np — 2 triangles in the corresponding Delaunay triangulation if in the Voronoi tessellation 
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there are n, interior- and nj, boundary polygons. Scanning electron microscopy is explained and 
several SEM pictures of cellulose nitrate are given. AVS is used to do the pore perimeter calcula- 
tion by connecting the modules read any image with image measure and image viewer, the latter 
two of which are in turn connected together. To calculate the pore area connect read any image 
with sketch, histogram, and image viewer. Then connect sketch with histogram, sketch with 
image viewer, and histogram with print field. A 3-d process is a homogeneous Poisson point 
process if the number of points in any region of volume V has a Poisson distribution with param- 
eter AV and the random variables corresponding to the number of points in disjoint regions are 
independent among one another. Among the characteristics calculated are, in our terminology, the 


ee. f thi a 
perimeter of polyhedron P. = )7j'*; e%, the surface area A. = Oj", O72, Ap”, and the volume 


is Ve = , Se V,°*). His numerical results given (see also, Jackson et al, 1999) are ng = 
5.203390, n& = 39.191489, n? = 26.127660, nf = 15.063830, P. = 17.358196b, A. = 5.910997b?, b® 
being the mean cell volume; the coordination number is of course 4. The distribution of n? shows 
the distinct Voronoi characteristic of containing exclusively of even numbers, whereas that of n& the 
equally distinct characteristic of being exclusively in multiples of three. This has been explained as 
odd number implying an impossible fractional number of edges, there being 3/2 edges as there are 
vertices. For the cross section data he found n¢ = 5.831 and P, = 3.563s, s” being the mean polygon 
area. In membranes, if the diameter of a particle is less than that of a pore then the particle will 
entrain the pore unless it should encounter other fouling particles, but if the particle is bigger than 
the pore then it will be retained. He assumes fouling particles to be spherical. Ascribed to each face 
is the diameter or size of its largest inscribed circle. Algorithm 6.2 is used for this purpose; d is the 
largest inscribed circle diameter sought for. 


Algorithm 6.2 Find the largest inscribed circle 


d<0 
for all combinations of three edges of the face 
find the corresponding triangle 
find the largest inscribed circle of this triangle 
if the centroid lies within the face then 
if the perimeter lies within the face then 
if the diameter > d then 
d+¢ the diameter 
endif 
endif 
endif 
endfor o 


In carrying out his simulation, the mean particle diameter dy is calculated from the mean pore 
diameter d, by d, = kd, where k is a constant. The corresponding standard deviation of particle 
sizes is op = 0y(d,/d,). The distribution of particle diameters is dependent on both k and another 
dimensionless parameter m = (op/dp)/(ov/dy). Then the central limit theorem gives the frequency 


distribution and thus the diameter of the i'” particle as d; = dp + op (ed ri 4/2) /Va/ 12), 


where 0 <r; <1 are random variables and gq > 12. Interaction between particles and membranes 
starts from Algorithm 6.3 which has been adapted from Jackson (ibid.). 


Algorithm 6.3 Interaction between particles and surface pores. 


while another particle exists do 
ait; 
pi <— particle ; 
find particle entry coordinates; 
c; + of the nearest Voronoi cell; 
fi < the nearest face of c;; 
if the face is fouled then 
for all the neighbouring faces f; 
if f; is not fouled then 
{f} = fi 
endif 
endfor 
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else 
{i} efi 
endif 
if {f} is not empty then 
j<0 
while particle travels and j < |{f}| do 
fy ds 
if d, > d; then 
if the centre of p; lies outside the largest circumscribed circle of f; then 
fi < the neigbouring face f; of the adjacent cell; 


endif 
if dp > 1.1d; 
particle cakes or blocks surface; 
else 
particle blinds surface face; 
endif 
else 
particle enters membrane; 
endif 
endwhile 
else 
particle cakes or blocks surface; 
endif 
endwhile D 


Backflushing the membrane may clean blocking- but not blinding particles. Once inside the 
membrane, each particle independently and randomly walks until it can get out of the membrane or 
can go no further and thereby necessarily fouls a cell. The interpretation shown here as Algorithm 
6.4 is an adaptation of that given by Jackson (1994). 


Algorithm 6.4 Interaction between each internal particle and the pores. 


for each internal particle p; 
{c} pores to enter; 
jel 
while j < |{c}| and p; has not entered a polyhedron do 
find pore to enter; 
if the polyhedron is fouled then 
particle enters the polyhedron; 
if there is no space in the polyhedron then 
particle blocks or blinds pore; 
endif 
endif 
icitl 
endwhile 
endfor 
On the cake, the particles drop, roll, or nest among one another. The cake formation is 
described as Algorithm 6.5 which closely follows the description which he gave. 


Algorithm 6.5 Particle cake formation, cf Jackson (1994) 


oD 


while next particle exists do 
pi < next particle; 
if route not directly leads to surface face then 
pi and p, collide; 
find rolling path; 
if p; does not roll off then 
p; and pe collide; 
if p; does not roll off then 
p; and pz collide; 
p;, comes to rest; 
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endif 
endif 
endif 
endwhile D 


Fifty particles is the empirical rule of thumb for the number of neighbouring particles to 
consider at the surface and cake level. Good care has been taken by his program to ensure that no 
two particles within this set overlap. All of his particles are gentle for they bounce not but often roll, 
for these no hard billiard balls but mathematical particles with assumptions. In this way a spherical 
particle would first vertically drops, then touches another particle and starts rolling off the latter 
along its surface. Once the two particles are abreast with each other they part. A dropping particle 
can either roll off one- or simultaneously two particles at a time. Rolling simultaneously off three 
particles is not physically possible. For if we picture i meets 1, rolls off 1, meets 2, starts to roll 
simultaneously off 1 and 2, and then meet 3. At this point it can either choose to give up 1 to roll 
simultaneously off only 2 and 3 or, it nests on the cradle of 1, 2 and 3, and thereby stops. There are 
three phases or layers associated with this membrane packing model. These are the cake-, surface- 
and internal layers. For the internal layer the thickness is limited by the membrane thickness, for 
the surface layer by the maximum diameter of the dropping particles, and for the cake layer by 
their amount and sizes. The flux decline across the membrane model is calculated by Algorithm 6.6. 
Here C,, is the particle concentration, Po initial pressure, Py final pressure, J, initial volumetric 
flux rate, that is to say, clean solvent flux for membrane, d diameter of the circular membrane, n 
the number of particles to be dropped on to the membrane in this simulation, L the depth of the 
packing retion, dy the mean spherical particle diameter, : the solvent viscosity, €, the free volume, 
that is to say, the ratio between the volume of voids and volume of bed, v initial velocity of the 
fluid, P,; pressure between the cake and the surface fouling layer, R; the pressure-flux relationship 
for the particles in contact with the membrane, R the pressure-flux relationship of the membrane, 


Algorithm 6.6 Flux decline across the membrane model. 


define At, C,,, Po; 
ve (4/1) Im/C; 
t< 0; 
while simulation should last do 
tett+At; 
n + £(At, t, Jm)3 
find the packing density, dimensions and number of particles of the cake; 
find the packing density, dimensions and number of blocking- and blinding particles 
of the surface fouling layer; 
find the packing density, dimensions and number of blocking- and blinding particles 
of the internal fouling layer; 
P. © the pressure drop across the cake calculated from v = (Po= Pr) a a 
P,¢+ P—-P, 
P, + the pressure drop for the surface layer; 
P; < the pressure drop for the internal fouling particles; 
Ri — Pari 
Ro+ Ri +R; 
Fy its - P;/ Ra; 


endwhile D 


Fluz is the rate of flow by weight. The number of particles per second is obtained by 
multiplying} the concentration by the volumetric rate. And if we divide this by the area of the 
membrane we get the number of particles per second per unit area. The random walk of particles 
during removal by backflushing is shown here as Algorithm 6.7. 


Algorithm 6.7 The random walk of particles during removal by backflushing 
remove cake and all surface blocking particles; 
while nezt particle exists and backflushing do 


backflushing < false; 
if not blocking then 


{ It is possibly a typing error in Jackson (1994) when he says dividing instead of multiplying here. 
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remain blinding; 
if there are more particles in cell then 
blocking < true; 
else 
no more particle; 
backflushing < true; 
endif 
else 
blocking < true; 
endif 
if blocking then 
do 
if cell inlet pore unfouled or fouled and the particle passes through adjacent pore then 
move into the next adjacent cell; 
else 
remain blocking; 
backflushing < true; 
endif 
while the particle is still within the membrane enddo 
particle has left membrane; 
backflushing < true; 
endif 
endwhile o 


Filters are sometimes made of non-woven materials, in which case they are modelled as tessel- 
lation by random straight lines. Wilkinson et al (1986) studies this type of filter by modelling it as 
intersecting random rods. The position and orientation of a rod is defined by one point, an angle 
0 < 6; 180° and the diameter of the fibre d;. The free area is ¢4 = (A — }> A,)/A, where the total 
fibre cross sectional area is }) As = D0, lidi — Yom Yon Oi did; Cosec O;;, where 1; is the length of the 
fibre, 6;; is 1 if lines i and j intersects and 0 otherwise, 0;; the angle of intersection of the fibres. 
Their model worked well for small fibre diameters compared to those of the particles. 

Algorithm 6.8 rewrites their method in an algorithmic form. Here 0 < rj; <1 is a random 
number with uniform distribution, d, the critical diameter, that of a circle which just fits the inside 
an irregular polygon. 


Algorithm 6.8 Non-woven fibre simulation, cf Wilkinson et al, 1986. 


for all particles 7 do 
dj —d+o(d), Tr — 1/2)/.fR/12; 
choose the particle position 0<a2;< X andO0O<y;<Y; 
for each layer of the filter 7 do 
dj < the largest circle which fits inside the polygon; 
if d, > d; then 
the particle is retained at this position; 
endif 
endfor 


endfor 
o 


The fibre here is in the form of fibre matt, which stacks one upon another in layers. The 
equation for the volume ratio suspension to inlet of a fibre matt, V,/Vo = exp(—&n) where € is 
the mean capture efficiency of the filter, is similar to that of the concentration ratio in deep bed 
filtration, C/Co = exp(—Ka). The mean capture efficiency is related to the layer efficiencies by 
€=1-[(1-&)(1—&)---(1— €,)]'””, and is dependent on d, a, €, ds, of, where the subscript 
f means filament. 

Jafferali (1995) studies Voronoi tessellation and applies it to microfiltration. The program 
which he used to create the Voronoi structure has the algorithm of Algorithm 6.9. 


Algorithm 6.9 Construction of 8-d Voronoi network (cf Jafferali, 1995) 


find set S which contains m nuclei nearest to the nucleus point p; 
find the nucleus q which is nearest to p; 
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find a nucleus r € S,r # q, such that the triangle pgr has the minimum circumradius, 
that is one facet of the Delaunay triangulation; 
do 
find s € S,s#r#q such that the circumradius of the tetrahedron pqrs is minimised, 
that is a Delaunay tetrahedron 
find m,n,o€ S such that m#4nZ4oAqHF¢r F58 and the tetrahedra pgrmj, pqsn 
and prso all have minimum circumradii, these are also Delaunay tetrahedra 
until each triplet pp;p; is in two Delaunay tetrahedra, that is every face of every 
Delaunay tetrahedron containing p as a vertex is closed enddo 


To find the centre of the circumscribed circle of the triangle ABC, let IL4, Ig and Ig be the 
planes perpendicular respectively to the sides opposite to A, B and C’. Respectively let n4, np and 
nc be the vectors normal to Il4, Ig and IIg, and a, b and c the position vectors of A, B and C. 
Then II4, Ig and Hcg will intersect at a point the position vector of whom is r which is obtained 
by solving the three simultaneous equations (r — n;)-a = 0 where i = A, B, C, ng = no X na, 
na =b-c, nc =b-a,d=0.5(a+b), e=0.5(b +c). It follows that rz = [m1 — (regay + 124z)|/az, 
ry = me + mer, and rz = (m4 —mg¢)/(m7 — ms) where my =r-a=a-ng, m2 =Tr-d=d-no, 
Mg =7r-e = €-Na4a, M4 = (Maz — Mdz)/(dyaz — Gydz), M5 = (Azdz — dzaz)/(dyaz — dydz), 
Mg = (M34z — M1 ex) /(€yax — Ay€x), and M7 = (Azeg — €z4z)/(€yAg — Ay€z)-. The centre p of the 
circumscribed sphere is obtained by solving the four simultaneous equations the first one of which 
is (py — Ay)? + (py — Ay)? + (pz — Az)? =r? while the other three have B, C and D in place of A. 
Solving these we get pz = (m2—m4)/(m3—m1), py = M1 £+Mae, pz = [2pe(Ca — dz) + 2py (Cy — dy) + 
ez — €1]/[2(d, — c,)| where m1 = [(az — be) (dz — cz) + (da — Cx) (bz — az)]/[(Cy — dy) (bz — az) + (by — 
dy) (dz —¢z)], M2 = [(e4 — eg) (bz — az) — (€1 — €2) (dz — €z)]/[2(y — dy) (bz — az) + 2(by — ay) (dz — €z)], 
m3 = [(a@z — dz)(bz — Cz) + (be — Cx)(dz — az)]/[(cy — by) (dz — az) + 2(dy — ay)(bz — cz)], m4 = 
[(e2 — es) (dz — az) — (e1 — 4) (bz — €z)]/[2(cy — by) (dz — az) + 2(dy — ay) (bz —cz)], €1 = a3 + ay +a, 
eg = b2 +b) +b2,e3 =e +c +c, and eq = d;, + di + dz. Statistics obtained from his simulation 
include ne, nf, n¥, n£, perimeter of polyhedron, area of face, and volume of polyhedron. The 
number of faces and the number of edges are found during face construction the method of which is 
Algorithm 6.10. 


Algorithm 6.10 Find faces of a polyhedron, cf Jafferali (1995). 


find the set of all the n vertices of the polyhedron, {v}; 
for all combinations of three vertices vj, vj; and vz of {vu} do 
find the plane II which contains v;, vj and vg; 
if all the remaining vertices lie on one side of II then 
Ui, Uj, Uz together with all the other vertices which lie on IT form a face; 
endif 
enddo D 


To find out whether all the vertices lie on one side relative to a plane, find whether the sign of 
the distance from each of these points to the plane > changes. This distance is d = pa-n = (a—p)-n 
where n is the unit normal vector, n = pa x be/|ba x be| = [(a — b) x (ce — b)]/\(a — b) x (ce — d)I, 
and a, b and c form a plane while p is the point in question. To arrange the vertices of a face in 
cyclic order, first pick one point among them and then find the angle made by the lines from it to 
every two combination of the remaining vertices, by 9 = cos—![(Up0j - UpD;)/|UpVi| - |UpT;|]. Two such 
lines which maximises the angle are both edges, and the angle from either one of them to each of 
the remaining vertices increases as we tread the edges from one vertex to another in succession. 

Nuclei points of the modified point process where d > dmin are generated by Algorithm 6.11. 


Algorithm 6.11 Nuclei points of the modified point process, Jafferali (1995). 


n; dmin; 

find pi; 

ae]; 

while i < n do 
tet=1; 
find p;; 
for j =2 to (i— 1) do 

find d(pi, p;)3 


if d(p;,D;) < dnin then 


i 


endif 
endfor 
endwhile 


1; 
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Dd 


The statistics that he found for these modified Voronoi structures are shown here again in 
Table 6.1. 


dmin 
0.0 
0.2 
0.4 
0.6 
0.8 
1.0 
1.2 


Ne 

226 
218 
215 
226 
229 
214 
204 


min(nf) 


8 
8 
8 
8 


9 
10 
10 


max(n¢) min(ny,) 


23 
24 


max(ny) 
42 
44 
44 
44 
38 
36 
36 


Table 6.1 Statistics of the modified Voronoi cells, cf Jafferali (1995). 
His asymmetric Voronoi tessellations have cell volume increasing with the z-coordinate position 
of the cell. The procedure for finding the dmin between two cells in such structure is described as 
Algorithm 6.12. Here the volume ratio is 6 = Vo/Vi, and thus 6 = yo/y +1, dmin = mz; + c where 
m = (y—y1)/(z — 21). 

Algorithm 6.12 Asymmetric Voronoi tessellation, cf Jafferali (1995). 


6; 


yl; 

yo + 63 

et 6/3; 
m<€(1—c); 


dmin <— k[(1 — ¢)z +e]; 


o 


The statistics that he found are given for reference as Tables 6.2 and 6.3. Real numbers have 
been rounded to leave at most four decimal points to make it easier to read (cf Jafferali, 1995). 


point distribution 


Poisson modified 
Ne 1276 1171 
do, n? 33340 29946 
ne 26.1285 + 6.4207 25.5730 + 3.5573 
,n§ 50072 44945 
ne 39.2414 + 9.6235 38.3817 + 5.3317 
v,ni 19284 17341 
ni 15.1129 + 3.2114 14.8087 + 1.7821 
Vi 0.60733 0.6582b° 
V. 4.76 x 10-408? +2 10-4 5.62 x 10-40° + 6.8 x 10-4 
>; pi 641.3870 664.1398b 
De 1.0053b + 0.2290 1.13436 + 0.1058 
>, 4: — 40.99300? 44.2189b? 
A. 0.03210? 0.037807 
Table 6.2 Comparative Voronoi statistics, Jafferali (1995). 

distribution 
Poisson modified 
min max min max 

V 88x 10750? =1.4x1079b® =3.9x 10-40? ~=8.4 x 10748 
ni 7 28 10 21 
n& 15 78 24 57 
ne 10 48 12 34 


c 


Table 6.3 More comparative Voronoi statistics, Jafferali (1995). 


asymmetric 
1069 
27749 
25.9579 + 3.5301 
41643 
38.9551 + 5.2740 
16032 
14.9972 + 1.7514 
0.591263 
5.53 x 10746? + 6.0 x 1074 
554.016 
1.0365b + 0.3255 
36.4954b? 
0.03410? 
asymmetrical 
min max 
9.9 x 10-58? ~=—3.7 x 107388 
10 21 
24 57 
12 34 


The Voronoi domain created has ragged boundaries due to polyhedra protruding. The method 
used in his thesis is to slice the domain by a horizontal plane and then redefine those vertices and 
edges on that plane and at the same time reject everything above it. This is more similar to a 
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carpenter filing away at a block of wood than a stone mason choosing his stones. The volume 
of voids was updated by using AVS which does this by counting the number of pixels contained 
within a given area. Statistics of an asymmetric Voronoi structure were given. His study concerns 
with the simulation of dead end filtration and cake formation. The particle size distribution is 
quantified in terms of the ratio of the mean particle diameter to the mean inlet pore diameter, 
a, and that of the standard deviation of the particle diameters to the mean particle diameter, /. 
The first parameter indicates the relative size between the particles and pores while the second one 
that amongst particles. Thus a > 1 would mean that the particles are bigger than the inlet pores, 
while 6 = 0 that particles are mono dispersed. The interaction between the particles and cake 
is essentially the same as that used earlier by Jackson (1994). The evaluation of pore properties 
described is equivalent to Algorithm 6.13. 


Algorithm 6.13 Pore property evaluation, cf Jafferali (1995). 


identify edge pores; 
identify inlet- and outlet pores; 
for all pores do 

find the largest inscribed sphere; 

for all faces of the pore do 

find the largest inscribed circle; 

endfor 

endfor D 


Identifying edge pores amounts to first finding neighbouring polyhedra, and then finding for all 
faces whether each of them belongs to one and only one polyhedron. Finding neighbouring polyhedra 
amounts to finding for all polyhedra pairs whether each of them possess no less than three common 
vertices. 

Particles pass through a face into a pore if they could, but since they are spherical whereas the 
pores and faces are polyhedral and polygonal, the largest inscibed circle of the facet is computed. If 
the particles are irregularly shaped, they may fit through a polygon when challenging from a certain 
direction but not another. However, the assumption of spherically shaped particles provides most 
authors on the subject with an approximation both satisfactory and convenient, because then there 
is no need whatever to rotate them. Jafferali (1995) applies an algorithm equivalent to our Algorithm 
6.14 to find the inscribed circle for all faces. Essentially this involves computing an inscribed circle 
of a triangle arisen from every one of the possible combinations of three sides of the face. Then the 
largest one of such circles which lies inside the face is the inscribed circle of the face. Here a, b and 
c are the positional vector of A, B and C respectively; n is a unit vector from A to B, q normal 
to the plane ABC, and m perpendicular to both; s; are the sides opposite (x;,y;). From 3-d each 
face is transformed into 2-d in such a manner that A = (21,41) = (0,0), B = (22, y2) = (0,y2) 
and C = (#3,y3). From c— a = 43m + y3n, x3 and y3 are obtained; z — for Zentrum — is the 
centre (xz, yz) of the face in 3-d or its equivalent transformed centre (#,, %,); € is the nucleus of our 
polyhedron while €; its neighbouring nuclei. Then the k* combination gives the largest inscribed 
circle the radius of which is r,z. 


Algorithm 6.14 Inscribed circle of faces, cf Jafferali (1995). 


for each face do 
for all the ?C3 combinations j of sides of the face do 
R+<0; 
find A, B and C, the vertices of the circle formed by them; 
n+¢ (b-a)/|b—al; 
m << solve m-n=0,m-N =0, |m| =1 and |n| = 1; 
q+ (b—a) x (c—a); 
{z1, Y1; x2} - 0; 
y2 + |b—als 
y3 © [(agny + y3ny)mz — (agmz + Ney3)My] /(NyMg — NeMy)} 
r3 € [(z3mz + nzy3) — y3Na] /Mz3 
pr ae $43 
Se p/2; 
4 1/2 
ry (ST a(S— sd] /S; 
£z © (8123)/p = (8123 + $241 + 83%2)/p; 
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Ez « (s1y3 + 83y2)/p = (s1y3 + S2y1 + $3Y2) /p; 
LZ, — £,Mz + EzNg3 
Y2 — £zMy + FzNy3 
2, + ABC(2xz,y2)5 
de |z— |; 
CH]; 
for all €; do 
di — |z — &il3 
if d; < d then 
C+ 0; 
endif 
endfor 
if ¢=1 then 
ifr > R then 
Rer;; 
ke jj 
endif 
endif 
endfor 
endfor o 


Algorithm 6.15 finds the inscribed sphere of each polyhedron. A tetrahedron is formed from 
each combination of four planes of polygonal faces. The largest of all the inscribed spheres of these 
tetrahedra is the inscribed sphere of the polyhedron. The intersection of bisecting planes Ni4_ Ili; 
is the facet in (3 — kmod(3)) dimensions which is equidistant from II; and II;, 7 = 1 to k. Finding 
Il = ax + by+cz +d =0 amounts to solving A-(B x C) for a, b, c and d. Finding intersection of 
planes amounts to solving their equations simultaneously. Here II; is the plane containing the i*” 
face, II;; the plane bisecting I; and II;. 

Algorithm 6.15 Inscribed sphere of a polyhedron, cf Jafferali (1995). 


for each polyhedron do 
R+¢ 0; 
for all the "C4 combinations j of the polyhedron do 
for i= 1to4 do 


Atv? —v}; 

Bev} - 0}; 

Cevi—vup; 

II;(a;, b;,c;, d;) < solve A- (Bx C) =0; 
endfor 


find Thy; i= 2, 3, 4; 

pao Ths 

ry = d(pi, Th); 

ifr; > R then 
R153 
ke jj 

endif 

endfor 
endfor D 


Particles move down vertically in discrete time steps towards the membrane. Upon reaching 
the latter, the i*” particle is greeted by an entrance the aperture of which is circular with radius 
re. If ry < re it penetrates into the pore, ifr. <r; < 1.1r, it blinds and if 1.1r. < r; blocks 
the entrance. Blinded particles are not removed by reversing the flow. They become a part of the 
membrane. Reaching the membrane is by no mean their destiny, and the particles only begin their 
long and tortuous pelerinage hereafter by embarking on a random walk in the direction towards the 
centre of the earth. Each particle enters a pore via the face furthest away from the latter, and it 
leaves via the face closest to it. The position of a facet is that of its inscribed circle. If the lowest 
facet is not viable, the particle first repositions itself precisely at the centre of the inscribed sphere, 
and then either moves out from some facet lower than itself or, when all possibilities of travelling 
having been exhausted, become a residence of that pore. The repositioning part above may seem 
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like a gross approximation, but there is no reasons why this should not make a sound assumption 
if we consider the fact that a real particle is never spherical in the first place. But if the particle 
blinds or blocks a facet, then the calculation used by Jafferali (ibid.) becomes more accurate. That 
isa =£+X(z—&), where A = |d|/(|d| + Jal), lal = (r2 +172)? and d= | — z| —|a|. To summarise, 
a particle is in a perpetual search for a lowest facet, which of course is closest to the centre of the 
earth. This kind of study is important because it lets us know which part of the membrane is prone 
fouling by the particles. 

Non-woven fibres are simply stacked layers of material. Extra holes produced by needles which 
are used to increase their strength are difficult to model and are normally assumed away before 
the simulation (cf Chan, 1990). Layers are considered as flat when the fibres which make them 
are considered to be flexible. Capturing of particles must take into account the separation between 
the layers as well as the usual aperture perpendicular to the flow direction. This is an example of 
man-made stacked 2-d layers. Examples in nature are numerous, including the honeycomb and cell 
growth in tissues. 

In his simulation, Chan (1990) represents random lines with (y — y;) = tan6;(a — x;), where 
0 < 6; < a, a method of interior randomness. Other possible methods include the w randomness 
where the line is represented by x sina; + ysina; = d;, where d; is the distance of the line from a 
pole point p; which is usually taken as the centre of the area, and surface randomness where the 
position on the boundary of the area is chosen. 

Flow channels within a medium follows Poiseuille’s equation for flow within capillaries, Q = 
APrr*/8uAz, where Az is the length of the capillary. Darcy’s law relates the flow of a fluid 
through a porous medium to the overall pressure drop, AP = pLug/K, where ug is the superficial 
fluid velocity, E the depth of the medium and K the permeability coefficient. The permeability 
coefficient for packed media, K = «/ [K'S?(1—)*], where < is the voidage of the medium, S' the 
specific surface and K’ the Kozeny constant. Particles can blind during both the filtration and 
the backflushing stages. The backflushing efficiency is the ratio between the volume of the particle 
flushed away and the total volume of particles captured in the media before backflushing. 


§ 6.1 Separation processes 


The largest aperture in the screen is the determining factor by which the selection is made, in 
other words it defines the size of the cut. This is a more logical argument than the one by which 
it is the mean aperture that defines the cut size (cf Rose and English, 1973). This largest aperture 
determines the smallest size of the particles retained on the screen. But in practice much smaller 
particles than these will remain on the screen for various reasons the most important one of which 
is blindings caused by particles with sizes dmin < d < 1.ldmin. The factor 1.1 was first used by Rose 
and English (ibid.), and arises from the angle of friction y being approximately 20°. 

On a sieve at each instant, the area that is free from blinding material is A = A, — G(mo —m) 
where A, is the sieve area, G the area blinded by a unit mass of blinding material, m and mo the 
mass of respectively the blinding- and initial blinding material on the sieve. Also, dm/dt = —kmA = 
—km [As — G(mo — m)] = —am(8 +m), where a = kG and 6 = [(A/G) — mo]. 

In surface straining filtration the particles are larger than the pore size, in depth straining 
filtration both the pore and the particle sizes are commensurate to each other, and in adsorptive 
filtration the particle size is smaller than the pore diameter. Adsorptive filters can be made to 
have higher filtration efficiency, higher capacity and higher flow rate than depth straining filters 
(Raistrick, 1986). 

Examples of solid-liquid separations normally found in Chemical Engineering are filtration, 
sedimentation, flocculation, centrifugation, electro-osmotic and electro phoretic dewatering and hy- 
drocyclonic separation. Poole and Doyle (1965) listed some of the aspects in filtration which they 
thought need further investigation: the effect of rapid pressure increases on the approach to equi- 
librium porosity in cakes, for instance in rotary filtrations; the relation between drag forces on 
particles and particle arrangements and shapes; migration of fines within cakes and media; the effect 
of changing flow paths during washing on porosity. 


§ 6.2 Dead-end filtration 


Because membranes are porous media, the flux across them follows Darcy’s Law, 7 = Ap/(uR), 
where j = dV/(Adt), R is the hydraulic resistance, R = Rm + R- where R,, and R, are membrane 
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and cake resistances respectively. Assume R,,, = 0, then R = R, = am/A, where a is the specific 
resistance of the cake, m the mass of the filter cake and A the septum area. 


§ 6.3 The centre of gravity 


The formula for the centre of gravity of objects in general is cy = f rw(x)dx/ f w(x)dx. For 
a real object the theoretical procedure is complicated, since there is always the possibility that the 
density is not uniform, but the practical procedure is simple and straight forward, that is by simply 
hanging the object by strings in various positions and then find the intersection between the lines 
of string, provided of course that this is possible with the object. 

The centre of gravity is important whenever there is a gravitational interaction with an object. 
Both the translational and rotational motions of objects through space are relative to this point. 
The c.g. of a triangle is positioned at one third its height whereas that of a half circular disk is 
4r/3x from the straight boundary line. A c.g. always lies on the lines of symmetry when these exist. 
Other names for c.g. include geocentre, centroid and barycentre. 

The centroid of a triangle lies at the intersection of its median. That of a tetrahedron lies at 
the intersection of all the lines joining the vertices with the centroids of their respective opposite 
faces. Its coordinates are the mean coordinates of the four vertices. 

To find the centroid of a polygon, first tessellate it into triangles and find the centroid of each 
one of them. Then we have another system of point masses located at the centroids of these triangles, 
with the mass proportional to their respective areas. 

The centroid of a rod is at its mid point, and the weight of a rod is proportional to its length. 
This can be helpful when we want to find the centroid of a network of rods or tubes. 

The centroid of a quadrilateral is the intersection between its two bimedians. This is also the 
mid point of the lines joining the mid points of the two diagonals. 

The medians of a triangle which has its vertices at (a1, b1), (a2,b2) and (a3, b3) are the lines 
connecting these three vertices to the mid points opposite to them, the first one to (c,,d,) where 
C1 = (ag +a3)/2 and d, = (b2 +63)/2, and similarly for the other two. Their intersection is obtained 
by solving their equations, (d, — 61)/(c, —a1) = (y—b1)/(# —a1), and so on. This gives the solution 
as being the average value of the vertices, [(}73 a:)/3, (03 bi) /3}- 

Having mentioned the centroid, it is natural to add the other two points related to it. On the 
Euler line also lie the orthocentre, where altitudes intersect, and the circumcentre, where perpen- 
dicular bisectors intersect. 

Let a;; be (a; — aj). Then for the orthocentre we need to solve any two of the three equations 
of the altitudes. For example, —a23/b23 = (y — b1)/(# — ai) and (y — be)/(@ — a2) = —ai3/b13 can be 
simultaneously solved to give x = [ai(a2bi2 + agb31) — (a2a3 + b21b31)b32| /(agbe1 + a1bs2 + a2bi3) 
and y = — [a32(a1a2 — agag + beb31) + agi (aiae23 + b1b23)] /(a3be1 + a1b32 + a2bi3). 

Similarly the circumcentre can be solved from two equations, for instance (y—bi)/(t—a1) = m1 
and (y = be) /(x = a2) = mez where my = —a23 /be3 and ms = —ai3/bi3, to give t= b32(a2a3 + 
b2b3) /(agb2 — azb3) and y = ag2(aza3 + b2b3)/(a2b3 — agzb2). 


Figure 6.1 shows a 2-d Voronoi tessellation with its 
nuclei and centroids. Ironically, the centroid repre- 
sents the position of a cell better than its nuclei, be- 
cause it positions itself in the most balanced situation 
relative to the cell, whereas the nuclei does so with 
regard to its neighbours. The centroids are noticeably 
more evenly distributed than our Poisson nuclei. 


Figure 6.1 Centroid of the Voronoi tessellation. 


Because the centre of gravity is a property which has more to do with the individual cells than 
with their neighbours, there may be real situations in nature where the c.g. of individual cells in a 
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network has a role to play. For example, it is usually assumed in a mathematical Voronoi tessellation 
that all the nuclei stay fixed and do not move. But in a real situation this would imply that there 
are some nuclei which lie very close to the walls. As no such nuclei would be stable, they more likely 
move away, and in that case the most reasonable prediction is that they move towards the centre of 
gravity of the cell. 

Because they are stochastic in nature, computational VT could in theory take up any shape 
however awkward. Jafferali (1995), for instance, imposes constraints on the algorithm which creates 
VT’s, namely in the course of the generation no new nuclei can assume a position too close to those of 
the existing ones. Thus there is an excluded volume within the network which develops throughout 
the generation process, and the Poisson point process does not over the whole space. Again, this 
excluded volume is relative to the neighbours of the cells rather than to themselves. 

I suggest that an evolution centred around the centroids of a network can replace such constrint 
regarding the minimum allowable distance between two nuclei. And that this would be a more 
natural for a network to do because the nature of the procedure which makes it look inwards to 
itself, as opposed to being totally governed by its environment. An attractive feature of the movement 
of the nuclei towards c.g.’s is that the cells in effect takes into account not only themselves but also 
their surrounding, because, if nothing else, it was the latter who shaped their appearance in the 
first place. Therefore the nuclei would start to grow from a purely random position according to 
a Poisson point process. Then it would move towards the centroid of the cell, and as it does so 
the definition of its boundaries is also changed. The centroid is a stable position, as can be seen in 
Figure 6.2. Similar to the processes of covering, C’(V) introduced in § 6.3 on page 172, and dual 
Voronoi, Y"(-) in § 6.3 on page 172, this centroid process can be recursive, that is G”(-). However, 
unlike the other two mentioned, this process very quickly becomes stable; Figure 6.2 (b) is much 
different from (a), whereas it is very similar to (c). The difference is in the size distribution as well 
as in the location of nuclei. The nonuniformity in the shape of a half-circle which can be seen in 
Figure 6.2 (a) is propagated to both Figure’s 6.2 (b) and (c), but it has become much fainter than 
in the original. 


Figure 6.2 Nuclei evolution towards a centroid; (a) a Voronoi tesssellation based on the Poisson 
point process, V, (b) the VT around its c.g.’s, that is G(V), and (c) the same process applied 
again to give the second order G?(V). 


The program in § A.22 finds the centroid of a Voronoi network and G”(-). 

There are three main steps in finding the centroid of a polygon. For each triangle of the 
triangulation of the polygon, first find the coordinates of its centroid, then its area. And then the 
centroid of the polygon is the weight average of the centroids of all the triangles. The first step 
involves finding the equation of two medians, which can be done by finding the mid point of an edge 
and then the slope to it from the opposite vertex. The second step can be carried out using the 
Heron’s formula which finds the area by means of edge lengths. 

Without loss of generality, let the polygon has n vertices at (a;,b;), where 7 increases from 1 
to n in the clockwise direction. Assuming that all triangles which triangulate our polygon share 
one vertex at (a1, 61), so that they all are Aj;j,, where i = 1, k = j + 1 and the three vertices are 
in the clockwise direction (a;,b;), (aj,b;) and (a,,b,). There are three medians in each triangle, 
corresponding to the three vertices, but we only need two in order to solve for the centroid. Any two 
of them will give the same result, but in order to fix the algorithm let us choose those medians going 
out from the vertices 7 and j. Then, since there is no loop in the procedure, we can save the space 
by describing it here in linear steps as ¢; + (a; + an)/2, dj + (0; + bn)/2, mi + (di — b;)/(ci — a4), 
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and the first median has the equation (y; — b;)/(x1 — a;) = m;. Likewise the second equation can be 
obtained as (y; —b;)/(x1 —a;) = m,. Here the subscript of the centroid coordinates tells us in which 
triangle it belongs. As the result of the above, we obtain for the k" triangle x, = (a; + aj + ax)/3 
and yz, = (0; + b; + by) /3. 

The second step is to find the area. For this, we first find the edge lengths s,1 ¢ (a2, + b3,)'”, 
8h2 (a5, +B)? and 8x3 + (az; + eye then the half-perimeter 5% < (S41 + $42 + $%3)/2, and 
then the area Ay < (84 (8x — Sh1)(Sx — 842)(Sk — Sp3))*/?. 

The third step solves for the polygonal centroid (#,y) as & + (0, @rAk)/ > >, Ak and y + 
Oo, yrAn)/ >>, An- Notice that the first step is parallel to the second- but not the third one. 
Moreover, finding c’s and d’s are parallel, as well as finding the various s,;’s. So we could vectorise 
these using parallel processing. But on Matlab all the matrix operations are already vectorised, so 
we only need to put all the components to be run in parallel in a single matrix and then do all the 
operations at once as a single operation on the matrix, provided that this is possible. 

The general algebraic solution to the above algorithmic procedure is obtainable by solving all 
the equations involved. This gives for any polygon the following theorem. 


Theorem 6.1. 
de gla + aj + On) (Gidjn + Oj dni + anbi;) 
— Se ee eee (24) vi 
3 d, (aadjr + a; bp: + ay bij) 
and 
_ CE +b; + bp) (asbjx + ajbgi + anbi;) (25); 


3 do, (Gadjr + aj bn; + andj) 

Proof. Assuming that we accept the centroid of a triangle to be at the mean coordinates among its 
vertices, and the centroid of massive point bodies is their average coordinates weighted by their mass, 
then when k = 1 it is obvious by looking at the Equation’s 6.24 and 6.25 that this theorem is true. 
Let us suppose that both these equations are valid for 7 triangles. We can present both of them in 
the form n/d where, necessarily, n = }>,,x;A; andd = )°., A;. If we now add another triangle to this 
cluster, the new triangle will also be represented as a massive point and it will add the terms xA to the 
numerator and A to the denominator, making n/d = (n+ xA)/(d+ A) = (Qo(y41) tiAa)/ Vga Ai 
Now, « for a triangle Ajj; is (aj + aj + a%)/3, and A calculated from Heron’s formula above is 
(agbjn + ajbng + axbi;)'/? /2, which yields Equation 6.24 and 6.25 for « + 1. This proves both 
equations, and thus Theorem 6.1, by induction. oD 

These general formulae, Equation’s 6.24 and 6.25, cover the triangle itself, when we consider 
the latter as a polygon, and the solution when & = 1 in this case reduces to (2, y) being simply the 
average coordinates of the three vertices, which is the same as that we have earlier mentioned. 

Earlier we have seen the evolution of the nuclei towards centroids in two dimensions. Let us 
look at the same thing in three dimensions. For this purpose, the program of § A.21, which is used 
in § 6.5 and explained by Algorithm 6.18 in page 177, is still inadequate, as its list of faces contains 
considerable amount of duplicates. This would have increased the resources required to do further 
work here. Therefore the program has been adjusted to that in § A.23 for the described in the 
following. 

Similar to the idea shown in Figure 6.2, we draw in Figure 6.3 the VT’s in three dimensions 
whose nuclei evolve towards the respective c.g.’s. 


(a) (b) (c) 
Figure 6.3 Evolution of nuclei towards centroids in three dimensions, namely (a) V?, (b) G(V?) 
and (c) G?(V3). 
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This shifting in the nuclei positions preserves the mean of the size while cutting down its 
variance. The pictures in Figure 6.3 show how a 3-d VT transforms into its first- and then second 
order adjustment. If we plot the cell size distribution, we will come up with a picture similar to 
Figure 6.4. The appearance of the distribution where sizes are discretised into bins like this depends 
on the number of bins chosen. 
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(a) (b) (c) 
Figure 6.4 Change in size distribution as a result of a nuclei adjustment process, (a) the size 


distribution of a V3, (b) G(V?) and (c) G?(V3). From a view point at 96° azimuth and 0° 
elevation. 


Figure 6.4 serves as a graphical presentation, but the statistics are the following. For our V°, 
G(V?) and G?(V3) respectively the cell sizes are 0.9345 +0.4025, 0.9495 + 0.2652 and 1.0315+ 0.2086; 
the third central moment 0.0332, 0.0059 and 0.0061; the fourth central moment 0.0936, 0.0137 
and 0.0036. These means are only meaningful when considered relatively with each other, since it 
depends on how one chooses the normalising basis. Here the basis chosen is the volume for the 
original Poisson point process, which is 1, divided by the original number of generators, which is 
400. When we built this V3, 236 cells had been disposed of as they were thought to be those along 
the boundaries; when G(V?) a further 107 and finally when we built the G?(V*) 46 more. This leaves 
us with 164, 57 and 11 cells remaining respectively in our V°, G(V) and G?(V%), all of which are 
shown in that order as Figure 6.4 (a), (b) and (c). Of course a great part of these are hidden behind 
others, so you can not possibly see them all there. 

The sizes in Figure 6.4 are based on a volume magnification factor of 1,000, which is corresponds 
to the increase in size by a factor of three, but those in Figure 6.4 are based on a magnification 400, 
which corresponds to a factor of (400)!/? increase in the size. 

As the new program has considerably changed from the old one, it is listed again in § A.23 
which has evolved from that which is listed in § A.22 on page 266. Wherever the sign © appears in 
this monograph it means copyleft not copyright. 


§ 6.4 Molecular dynamics 


The Lennard-Jones function can be written u(r) = —A/r® + B/r!”, where A and B determine 
respectively the attractive and repulsive parts. Let the range parameter o = (B/A)'/® and the 
energy parameter « = A?/4B, and the equation becomes u(r) = 4e [(a/r)!” — (o/r)§]. In other 
words, o is the diameter of one of the atoms and ¢ is the well depth, i.e. energy constant. The force 
here is not F = —du/dt = (24e/r?) [2(a/r)'* — (o/r)®], but F = —du/dr. 

Another approach in molecular dynamics simulation is to use numerical methods on the New- 
ton’s equation of motion. The most widely used is perhaps the Verlet algorithm, shown here as 
Algorithm 6.16. The initial values are rg and 11. 


Algorithm 6.16 Verlet algorithm, cf L. Verlet (1967) 


for i= 1 ton do 
find fj; 
Tiga 274 —Ti-1 + fiAt/m + O(A?*) and 
vi © (rina — Ti_1)/(2At) + O(AP). 
endfor o 
There are variants and modifications of Algorithm 6.16, for example the Leapfrog Verlet 
algorithm which has three steps instead of two, that is Up41/2 = Un—1/2 + fn/mAt + O(At*), 
Tri =Tn +Un41/2dt + O(At*) and vn = (Un41/2 +Un—1/2)/2 + O(At?); or the velocity Verlet algo- 
rithm where rp41 = TntUn Att frAt? /(2m)+O(At?) and Ungi = Un tAt(fntitfn)/(2m)+O(At*), 
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or in the form normally used in practice in which vj41/2 = Un + frAt/(2m), Trya = Tn + Un41/2At 
and Ung¢1 = Un + fry At/(2m). 


The van der Waals force is an attractive force acting between molecules. In the case of gases, 
each molecule consumes some space, so the dynamic volume is less than overall volume by an amount 
bn when 6 is a constant and n the number of molecules. Here b is Nav, where v is the volume slightly 
larger than the volume of each molecule of gas and N4 the Avogadro number, N4 = 6.022 x 10-28 
mol. The pressure of gas we see is the pressure of gas detected, which is less than the real pressure 
by F/A, where F = )° Fj and F; = >> f;. Here i runs from 1 to n and j from 1 to (n — 1) in the 
same set of molecules. The reduction of the total force comes from all molecules, but this reduction 
from each molecule is in turn affected by those molecules around it. Therefore i 4 j and j runs from 
1 to (n — 1) as mentioned above. Suppose that each molecule attracts another molecule by a force 
a. Then, in a given volume V, F = 7 Fi, = 0,00, fig = Ulm — Da/V = n(n — 1)a?/V?. For very 
large n we can say that (n — 1) & n, and therefore (P + an?/V)(V — bn) = nRT. 


The nature of this mutual attraction between molecules, which reduces the pressure in gas, 
comes to light in the case of solids. The distribution of charges around a neutral atom in solid 
fluctuates in the time scale r < 1071® s. This charge imbalance makes each atom behave as an 
electric dipole. This electric dipole has an electric field around it, which affects other atoms nearby 
and binds the two together. 


The electric dipole moment is p = qa, where a is a vector from the negative to the positive 
charge. The electric field around a charge q; has a magnitude E = q/(4meor”). It acts on a nearby 
charge q with a force of magnitude F = u = q.E = qiq2/(4neor”), where u is the potential energy 
of the two charges. This force itself is the Coulomb force, and the electric potential around q, is 
u/qa, that is V = qi /(47eor). 


For an electric dipole, V = (1/(47€0))(q/r1 — ¢/r2) = q(r2 —11)/(4ne0rir2). Since a <r, a 
being of the order of the atomic size, rir2 ¥ r?. Let 6 be the angle between a and r. Then, also since 
a <r, we have ro —11; & acos6. Therefore we now have V & gacos6/(4meqr?) = pcos0/4reqr? 
and consequently E, = —OV/dr = 2pcos6/(4reqr?) and Ey = —OV/(rd6) = psin6/(47eor?). 


We neglect Ey for long-range interaction, and simplify E, to E, = a/r?. We shall hereafter 
use the terms dipole and atom interchangeably. If this electric field is produced by atom py, it could 
induce in another atom po = a£&,, where a is the molecular polarisability of the second atom. This 
second atom will have an energy of interaction with the first one, u = —p.E, = —aE? = —aa?/r®. 
This force is attractive. 


The repulsive force is more complicated, and is generally thought to arise from the Pauli 
exclusion principle and the coulombic repulsion of the electrons in the outer orbit. It is generally 
assumed to be either u = A/r!? or u © aexp(—r/p), where p is a range parameter (cf de Podesta, 
1996). 


The Lennard-Jones potential is often written as u = —4e [(a/r)® — (o/r)'*], where o is a range 
parameter that indicates the approximate size of an atom and € an energy parameter that indicates 
the strength of the interaction between atoms. In other words, it is the minimum value of u. 


2 


The Lennard-Jones potential between two 
atoms has the minimum value u=e at r= 
1.12250. Figure 6.5 is a plot y = —4(1/x2° — 
1/z'7), where y = u/e and x = r/o. The 


force which the potential curve of Figure 6.5 ost J 
acts on another particle is shown in Figure 8 

6.6. The force produced by a dipole is calcu- of 

lated from F = —du/dr = 24(e/c)[(o/r)’ — S| 
2(c/r)!?]. But in our units y = u/e and -0.5+ i 


x = r/o above the equation becomes y = 
24(1/x* — 2/213), which is plotted in Figure -1f 5 | 
6.6. 


Figure 6.5 Lennard-Jones potential. 
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Figure 6.6 shows that the force is zero at 
the distance « = 1.12250 away from the 
particle. Closer than this point the repul- 
sive force increases rapidly towards infinity. 
Further away from this point, however, the 
attractive force increases to a maximum and 
then gradually dies down. This maximum 
force occurs when dz/dy = 0, that is at 
aw = 1.24450. 

| Figure 6.6 Force corresponding to the Len- 


39 1 11 1.2 a 1.4 1.5 1.6 nard-Jones potential. 


In gases, the molecules travel past one another so fast that they hardly notice the wells of 
negative energy surrounding other molecules, let alone stop and rest there. But even here the 
attractive force produced by these wells is probably what give rise to the term an?/V? in the van 
der Waals equation. Solid molecules, in contrast, hardly have kinetic energy and therefore prefer to 
sit in such wells of their neighbours, which is the reason why solids hold and neither flow as liquid 
does nor disperse like gas. 

The most cohesive structures of solids are crystals, where the Lennard-Jones potential cul- 
minates in a minimum cohesive energy of the lattice. This cohesive energy is U = Nau;/2. Here 
u,; is the summation over all pair potential energies of Lennard-Jones type, uj = Be u(raj). Let 
ij = QiyTo, where ro is the nearest neighbour distance. Then U = —2eN,4 [(a/r0) Ae — (0/10) Az], 
with the lattice sums being Ag = )7(1/a§;) and Ai2 = )>(1/aj?). These lattice sums are calculated 
from the network structure in terms of infinite series that converge very quickly. 

To find the optimum value of a/ro we shall let 6 = o/ro and write differentiate dU/db = 
—2eN 4(6b°Ag — 12)" Ayo) = 0. As a result, ro = o(2A12/Ag)'/®, the cohesive energy per mole 
U = —(A2/2Aj2)Nae and the cohesive energy per molecule u = U/N4 = —(A2/2Aj2)e. 


§ 6.5 Problem definition and algorithms 


Walls in a 3-d Voronoi tessellation are isomorphic to bonds that link between its cells. Therefore 
we can redefine the problem of a particle passing through the hole in a wall, into the cell chamber 
and out through a hole in another wall, and so on, as that of the passage of a mathematical particle 
through a tube that links between two cells, into a cell and out through one of the other tubes. In 
reality the solid parts which make up the partitions have thickness, therefore the holes through the 
walls will have a nonzero thickness and the volume of the cell chamber is smaller than that of the 
otherwise mathematical Voronoi cell. 

Each vertex of a Voronoi polygon has three walls of the latter attached to it. The centre 
of gravityt of these three faces is calculated and linked together. This truncates the coigns and 
produces for each polygon its dual self. The next step is to link together the midpoints of all those 
edges of the polygon that have a vertex in common. The edges of this last polygon bound and define 
the void volume of the cell chamber. 

Next, the cross section of the hole within each wall is taken to be the area on the original wall 
which is bound on all sides by the second order covering lattice of that wall. 

The filtering membrane in this case is assumed to be isomorphic and homogeneous. The size 
of the particles is taken to be reasonably smaller than the void in general, so that the attraction 
between particles and the adsorption to the wall play a more prominent role than the physical 
blockages by individual particles. The particles are all assumed to be of the same size. 

Algorithm 6.17 is an algorithm to do filtration that is being developed. Here N contains the 
neighbourhood information 


Algorithm 6.17 Filtration in Voronoi tessellation 


find (vg, Ca) <— Voronoi tessellation; 


{ aka centre of mass, centroid. 
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find d; «+ Delaunay triangulation; 
find N. + cell neighbours; 
¢ + Cq which lie completely within [0, 1]°; 
v + VU, which belong to some c; 
im} 


My work on filtering membranes results in another Voronoi algorithm and program which is 
different from the programs listed in §’s A.3 and A.4. Because it is also written anew from scratch, 
this new algorithm has nothing to do with the previous two as regarding the data structure and the 
logic of its method. It is given in § A.21. 

The current convention for variables is this. A single alphabet or entity means a list, for 
example c is the cell list and va is the list of all vertices original created. Another example is }, a 
structure for bonds which contains the list of the cells of bonds, the number of vertices of the face 
represented by each bond, the list and the map of these vertices, the ordered list of the vertices of 
each face and another list similar to this but cyclic, and the number of cells connected to each bond. 
Similar to the structure of b is that of bdr, or bg in Algorithm 6.18, which is a list of the bonds 
along the border, each of which is connected to only one cell. 

Two entities xy makes x of y, for instance vc is vertices of cells, which is a structure that 
contains the number of vertices of each cell and the list and the map of all these vertices. I try to 
preserve the space by creating a short but easy to understand naming convention. Also, the lines 
are put together, which means that the number of lines of the codes is more than what appears here, 
some of the lines having six logical lines or more to them. But the structure of the program still 
remains intact, therefore it should be possible without much difficulty to compare the program in § 
A.21 with Algorithm 6.18 which explains it. 

Three entities xy? or y? are y x y matrix of x, where the latter relates two entities of the 
former, for example bcc is the bonds mapped on to a cell matrix. In the algorithm, v, and cy, are 
respectively the vin and cin in the program, the vn and cn there being verbosely described in the 
algorithm as the number of v and c. 


Algorithm 6.18 Voronoi data structure for the study of membrane filters. 


(Ug, v2) + find a 3-d Voronoi tessellation; 
Un + index of 0 < ug <1; 
Cn + index of c all the v’s of which are in vy}; 
U © Ua E Un} 
Uc — Us E Cn} 
t, «+ find a 3-d Delaunay tessellation; 
for all ¢; € ¢, which are connected to two cells do 
tH ti 
else 
ba & tis 
end 
becet; 
be ba; 
find Delaunay triangulation in (3 — 1) dimensions for all faces; 
order the vertices of each face; D 


§ 6.6 Simplified algorithm for filtration 


The first algorithm we shall now develop is a reasonably simplified one. But it turns out to my 
surprise that this simple algorithm links us to the continuum percolation of hard spheres or even of 
particles with irregular shapes. This is because of the assumptions that we shall make. 

One such assumption is that the variance of the filter cell size distribution is small. For this, 
Jafferali (1995) would probably have applied his favourite constraint code on a Voronoi tessellation. 
The approach which I developed uses the operator G(-) on a V (cf § A.23). It is more logical and 
therefore is what we shall use for now. 

The second assumption is that the particle size is reasonably smaller than the size of the smallest 
void of the original system. This is the assumption which seems to link us to the percolation of 
spheres in continuum, the space considered being that of each void. 
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Also, we assume that the particles are attracted towards one another and towards the walls 
by various interactions which may include the van der Waals force and the force from electrostatic 
interactions. 

Let @ be the angle that the line from the c.g. of a each cell to the mid point of a face makes 
with the horizontal plane. Then another possible assumption is that all particles prefer a trajectory 
which goes through a face which has the maximum 6. 

The third assumption is that the volume of each membrane pore is 80 per cent the volume 
of the corresponding Voronoi pore. The fourth assumption, blockages due to the clustering caused 
by the attrition forces between particles can occur within pores. This excludes cake formation and 
blockage at entrances to, or exits from the pore. 

The last assumption above is similar to assuming that clustering due to attrition can occur 
in stagnation regions. Since the flow of the suspension reaching the membrane is strong, i.e. the 
pressure high, there can be no clustering there. The flow from one pore into another is also faster 
than the flow within one. This is amount to assuming that the effective diameter of the hole within 
each face of a pore is considerably smaller than its diameter. 

Upon reaching the membrane, each particle in the troop seeks out the bond closest to it, and 
begins its journey through the membrane. Within the membrane, it always follow the bond which 
has the greatest gradient available. Therefore the top bonds should map to the bottom ones one to 
one. But because the possible blockages during the course of operation, this may not be so and the 
mapping is instead one to many. 

First we consider the case in two dimension, so volumes becomes areas and the volumetric flow 
rate the distance travelled. In our discrete time, if v is the flow velocity, n the number of particles, 
py the total volume ratio of the particles and r their radius, then we have nmr? = vAt. Assuming 
that v = 1 ms~! and r = 100 pm, then n = 10°, /7 or approximately 318, particles for each time 
step. 

Next we will study the suspension in a square box in order to find out py. Because we shall 
assume that the particles flowing through the pores percolate as though there is no flow, we will 
consider here the suspension that is simply contained within our box without moving about. I feel 
that the assumption that particles could percolate that way in pores is justified since the flow is in 
steady state, protected from the turbulence outside by all the solid structures making up the walls 
of the membrane. 

The percolation program first generates random positions of the particles, then moves apart 
those which are too close together so that in the end they only touch each other. Particles which are 
separated by a distance less than 1.4r move towards each other until they touch. Lastly, touching 
particles never separate. 

We know that the total volume ratio of the particles has the upper limit of 0.9069, because 
that is the density of the closest packing of circles on a plane. 

Packing circles on the plane becomes densest when the circles are arranged as a hexagonal 
lattice with the packing density of 7/2V3. Packing of spheres is similarly at its highest density if 
the arrangement is that of the face-centred cubic lattice with the packing density 7/32. 


§ 6.7 Filtering problem when physical blockage is prominent 


Assuming particles to be spherical, the size of particles to be constant, and that this size is compatible 
with the size of the holes in the walls and the voids so that physical blocking is responsible for most 
of the blockages in the membrane. As in § 6.7, the hole in each wall is taken to be the polygon which 
results from recursively finding a covering polygon for the face twice. 

Redefine the problem of particles’ passage through walls and voids as that of particles travelling 
along edges of the dual lattice of the Voronoi structure, i.e. the Delaunay triangulation. Each of 
these edges corresponds to a face in the original physical lattice. To each edge is thus ascribed the 
details of the cross section of the hole through that face, namely the shape and the dimension, as 
well as the gradient it makes with the horizontal plane. Upon reaching a vertex, the particle ball 
will choose the next path, i.e. bond in the dual lattice, which has the maximum gradient to pass 
through. However, if this bond is too small or if it is blocked, the particle with choose from among 
the remaining paths the one which has the maximum gradient, and so forth. If the size of the particle 
is approximately that of the hole it tries to pass through, within two per cent of the latter, say, then 
the particle will blind the passage at that point. The difference between blinding and blocking is in 
the degree of tightness that the particle sits in the hole, which reflects in the degree of difficulty to 
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remove it by backflushing. This degree is not constant but a function of the relative size between 
the two parties involved, i.e. the particle and the hole. 


§ 6.8 Percolative filtering with very small particles 


At this point a new filtering algorithm is considered and investigated. We introduce an assumption 
that the particles are all of the same size which is very small compared with the size of the pores. 
So there is neither cake formation nor blinding by a single particle. Assume that the solid particles 
suspended in a fluid medium, being dragged downwards under their own weight. 

Assume that the van der Waals force plays a significant part and, similarly to § 6.8 on page 179, 
that these interactions between particles are governed by the Lennard-Jones equation. Assuming 
that each particle is spherical and acts as a dipole with the electric dipole p = 10e, where e is the 
electric charge on a proton, e = 1.602 x 107° C. Then (ef de Podesta, 1996) 


Pp 


7 Ager (26) vi 
and the energy of interaction between two particles becomes 
2 
ap 
u=—ak? = ——_, 27) vi 
” ~ Gnyenrs ans 


where q@ is the molecular polarisability of one particle under the influence of the other. If we assume 
that the repulsive force acts in such a way that the repulsive energy is u = c/r'*, where c is a 
constant, the Lennard-Jones potential becomes Equation 6.28. 
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Comparing Equation 6.28 to the form 
A B 
Up = ae + apy, (29), 


we have A = ap?/((47)*e2) and B = c. Or if we compare it to the form 
a G6 _ (Fy 12 ; 
uy = —4e [(2)° - (2) J, (30) 


then we have o = (B/A)!/° = (16cm?e2/(ap’))'/® and ¢ = (A?/4B) = a?p*/(4c(47)*e§). Notice 
that here o and ¢« are simply parameters, not the Stephan-Boltmann constant and the dielectric 
constant, o is arange parameter and approximates the size of an atom while ¢ is an energy parameter 
which shows the strength of the interaction between particles. At r = o the value of u, is zero, 
whereas u, has the minimum value of —e. 
The force between two particles is Equation 6.31. 
du Gap? 12c 
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At the equilibrium separation, ro, f = 0 and therefore Equation 6.31 yields the minimum 
distance in Equation 6.32, 
1 
32cm765 \ ® 
ro = ( B) 2) : (32) i 


The equation for the minimum energy is obtained by substituting ro from Equation 6.32 into Equa- 
tion 6.28, which gives 
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Returning to Equation 6.31, the attractive force has the maximum value at the point where 
dF/dr = 0. Differentiating F in Equation 6.31 with respect to r, we arrive at Equation 6.37. 


dF A2ap” 156c 


<— = —<F _- 87)vi 
dr (4m)2e2?r® or !4 sy) 
From this the distance where this maximum force occurs is given by Equation 6.38, 
tp 156e(47)*e5 _ A16cn*e9_ (38),; 
mm 42ap? Tap? ae 
Then by putting Equation 6.38 into Equation 6.31 we have this maximum force in Equation 6.42. 
6ap? 12c¢ 
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Let each of the spherical particles has a ra- 


mg dius r = 1 pm and the density p = 3,000 

Y kg -m7~* (compare the density of silicon, 

(a) which is 2,329 kg- m3). Then the mass 
& Tin 2 of the particle is m = pV = 3000 x (4/3) x 


(10-6)? = 1.26 x 10-™ kg and the weight 
is mg = 1.26 x 10-“ x 9.8 = 1.23 x 10-8 


N. We may now find the attractive dis- 

¢ i tance, the maximum distance whereby two 

particles will come together under the van 

i der Waals force. Figure 6.7 shows the cap- 
, Tm turing of one particle by another when one 


particle is fixed in space and another par- 
ticle has zero velocity but is free to move. 
Recall that the capture occurs when the ac- 
celaration due to the weight of the particle 
(b) mg g equals the acceleration a due to f. 


Y Figure 6.7 Force balance with one particle 
fixed. 


If we suppose that our particle have the same polarisability as that of a benzene in its gaseous 
state, CeHg, i.e. a = 11.61 x 10~*°, then in this case a = 11.61 x 10-*° F-!m* (ef de Podesta, 
1996) and Equation 6.42 gives us the separating distance between the two particles which gives the 
maximum attractive van der Waals force, which is 444 wm. Giving our particle other values of a, 
with the value of a for methanol gas CH3OH, i.e. a = 3.860 x 10-49 F-!m?*, we have this distance 
r = 379 pm, and with a = 1.647 x 10-40 F-!m’, that of water vapour, the distance becomes r = 
336 pum. 

These values are rather large compared with the radius of the particles, therefore we shall opt 
instead to a bigger size of particles, when r = 5 wm. With the radius of five microns, the particle 
has a mass of 1.571 x 10-1? kg and its weight becomes 1.64 x 10-'' N. Then if we adopt the value 
of a for water vapour, a = 1.647 x 10-49 F-!m*, then from Equation 6.42 the distance where the 
force is maximum becomes r = 167 pm. 

But this is not everything. So far we have only considered what two particles will do when 
one of them is fixed and the other one has no initial velocity. There are two other things that can 
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happen, the particles may both be moving down beside each other under gravity and the path of 
the particle which is effected by their mutual attraction. Even when two particles do not come 
together, their path can be deviated by the van der Waals force. But here for simplicity we shall 
neglect this and assume that particles either come and stay together or they experience no mutual 
force whatever. Only if they come together will their final path and velocity be affected, and these 
depend on their combined momentum. 

When particles move relative to each other, their captive velocity depends on their relative 
velocity. They tend to join each other more easily if their paths are along side each other and goes 
in the same direction. In the extreme case where both particles move in the same direction with 
zero relative velocity, they would come together across a vast gap in between indeed, had there been 
no frictional force due to viscosity that acts to drag them. 

In the present study, the friction due to the fluid in the medium is neglected, together with 
the relative velocity between the particles, for the purpose of deciding whether particles will come 
together. It is expected that on average particles coming to within the capturing radius of each other 
will come together. This radius is larger than the radius of maximum force in Equation 6.42 above. 
In other words, we expect our particles to behave like solid spheres with a well defined boundary for 
their sphere of influence and a much simplified capturing mechanism. 

We shall define this radius of captivity, r., to be one third that of the radius of maximum 
force rm. From Figure 6.5 this is the point where r, = 1.50, and because r,, is here 1.2445 we have 
Tc/Tm & 1.2. For Figure 6.5, this is the same as saying that r./ro = 1.34. 

One additional point is that, in real situation when particles form a cluster their collective 
value of the molecular polarisability will change, and this value will not be the same for the cluster 
as it is for the individual particles. But here we assume that they are the same, which means that 
the capturing radius of a cluster will be the same as that of the particles which form it. 

It may worth mentioning here that the molecular polarisability is related to an optical property, 
viz. the refractive index n,, by the relationship a = €9(n? — 1)/n, where n is the number density of 
molecules. For gases n = P/(kgT), whereas for liquids n = Nap/m, m being the mass. 

Next we shall concern ourselves with clusters of particles thus formed. All clusters will be 
assumed to be a closest-conglomerate of particles which form them, with the densest packing density 
possible in three dimensions, that is n/3/2 of the face-centred cubic lattice. Numerically this is 
0.7405. 

Notice that the packing density in two dimensions can be higher than this. Packing circles on 
a plane is the densest of all with its packing density of 7/ 2/3. 


2r 


Shown in Figure 6.8 are the hexagonal lat- Nr ae 
tice packing in two dimensions and the cu- [ \ 
bic close pack. In the cubic close pack- ie A / 
ing spheres in every third layer lie verti- . 
cally straight on top of one another. Each 
face of its cubic section looks like the sec- \ J 
ond picture in Figure 6.8. Similar to the re Pe S 
cubic close packing is the hexagonal close ( | | 
packing spheres in every alternate layer of \ J F ae 
which lie over one another. Both the cubic- (a) Se = 
and the hexagonal close packing have the a 
same packing density which, in the case 
of the former, is calculated, from the sec- a4 } 
ond picture, as p = >> vs/ >> vc, where vs Se wf reel tae 
are the total volume of sphere segments 20% ( 
in the unit cell which has the volume v,. Pa p> 
Here uv, = (2V2r) and ov, = (8(1/8) + [ ant 
6(1/2))4r?/3. In 2 dimensions, the dens- ] 
est packing is calculated from the first pic- es 
ture. (b) 


12 
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Figure 6.8 The packing density calcula- 
tion of the closest-packed densities. 
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The density of the packing is }> a-/a,, where a, are the areas of the circles and a, is the area 
of a rectangle. These two areas are namely a, = 2(rr?) and a, = 2r(2V3r), which give the density 
p=1/2V3. 

If the particle size is very small there will be an effect from quantum mechanics. The de Broglie 
wavelength, \ = h/p, exists for all particles large and small. Here p = mv, m = pV and V = xd?/6 
for a spherical particle. Figure 6.9 shows a graph between diameter and wavelength of particles. 


We can see that the quantum-mechanical ef- 
1? bs chan bans cucceseed fect when the pore size is in the order of mi- 
. a es ee crometre will become prominent if the size 
of particles challenging it is in the order of 
nanometre. When this is the case there will 
eras coe sgeeo 4 be not only the interference effect of the 
re ne eer te cer 1 particle-wave passing through slids but also 
ee : oS ee the effect of confined particles which pro- 
duces standing waves of these wave-particles 


10 Be es : 4 “cosa in each pore. These standing waves of the 
wel ees ate oi. particle-wave have A, = 2D/n. 
wo Figure 6.9 Quantum-mechanical effect in 
10° 10° 10” 10° r 
dim filtration. 


§ 6.9 Percolation within percolation 


Particles suspended in a fluid generally have their diameter dp much smaller than the pore 
diameter d, of a filtering membrane. Yet these small particles can cause blockages of the membranes 
due to attrition among themselves. Since in this case d, < d,, blockages due to blocking or blinding 
of the individual particles (cf Jackson, 1994; Jafferali, 1995) become out of question. 

Having investigated both the percolations of networks and continuum, I suggest that the block- 
age of these smaller particles in membranes is due to a double percolation phenomena, one the per- 
colation of the suspension continuum, the other the percolation of the centroidal Voronoi network. 
As a reminder of a centroidal Voronoi network, it is a Voronoi tessellation on generator points which 
are the centroids of a Voronoi network which either is generated from Poisson point generators or is 
another centroidal Voronoi network. 

Because percolation is a study of the behaviour of two phases, and because in general p, # 1/2, 
there are no less than three regions of behaviours to consider in each percolational investigation 
instead of two (cf Tiyapan, 1997, KNTs(ix); also in Tiyapan, 2003, KNT<a(iii)). When p, < 0.5, these 
three regions are p < pe, De <p < (1— pce) and p> (1—p-). When p, > 0.5, they are p < (1 — pe), 
(l—p.) <p<p,- and p> pe. The case where p, = 0.5 is assumed to be very rare in nature, and so 
can be neglected in the present study. 

When a suspension becomes so high that the average interparticle distance has become such 
that the attrition due to van der Waals force is prominent, the suspension will solidify into a moisted 
bed of particles. According to the percolation theory, we may define the point where this spontaneous 
solidification occurs to be that point where there is a single cluster, under a mutual van der Waals 
force, which traverses the whole continuum in a certain well-conditioned direction, that is to say, 
a direction which may represent the diameter of the network. Furthermore, let us call a critical 
concentration p, the minimum concentration at which this infinite cluster appears. 

Then we have for our suspension a continuum percolation with three regions of behaviour 
similar to those we have found in the case of network percolation. Furthermore we map the space 
of pe on to that of 0 < p. < 1, where p, = 0 means there are no particles suspended in the fluid, 
in other word p = 0, and p, = 1 is where the suspended particles form a bed in the closest packed 
structure, that is p = pmax. We also assume, without the loss of generality, that p, > 0.5. Then we 
have the following as the three regions, p < (1—p,), (l—p.-) <p<p, and p> pe. 

We are interested in neither the cases p < (1 — p,) nor p > pe, since the former implies too 
dilute a concentration for the attrition due to the van der Waals force to cause an infinite cluster, 
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while the latter means that the suspension is so concentrated that they solidify instantaneously, 
simultaneously in all pores. 

When (1 — p.) < p < pe, all pores has an equal probability of being solidified, and so the per 
cent total solidified cells now depends on the topology of the network, in this case Voronoi, and the 
probability where the critical phase change occurs becomes the critical probability of the network. 

In the case of dead-end filtration the flow is in one direction, therefore the critical flux reduction 
which is the result of percolation of the network occurs at the point where the cross section of the 
network, not the network itself, percolates. This is because such percolation in the cross section 
in the plane perpendicular to the flow direction will cause a bottle neck in the flow and therefore 
determines the flux. This phenomenon is summarised in Figure 6.10. 

Next we must define the critical probability of the overall system. Since the system involves two 
kinds of probability, that is continuum and network, and assuming p,, and p;, are two independent 
probabilities, then we have the overall probability is 


Pe = PoPea; (43)vi 
where (1 — pe.) < po < Peg- 


SB > 
0 ! Pp 


continuum 3-d 2—d section 


Figure 6.10 Percolation within percolation; 0 < po < pmax, 0 < po <1, Pe = Peg * Po- 
If our membrane is homogeneous, the reduction in the area perpendicular to the flow becomes 
AA = Ayp-,, where V, /Vi = Ay /Az, Vy and V; are respectively the void volume and the total volume 
of the membrane, and similarly for the areas A, and A. 


§ 6.10 The first part, suspended particles 


Because of the complex nature of the problem, there is no single algorithm but rather there is 
an algorithm for each job. The first task is to study the percolation of the spherical particles under 
the van der Waals force. 

Continuing from the development in § 6.10, particles are spherical in shape with rp = 5 ym and 
the capturing radius is r, = 167 x 1.2 = 200 wm. First we shall study particles within a cubic box 
of side length 2 mm. Particles start as a suspension with no obvious velocity. They stick together 
and to the walls. 

When particles come together, they form a porous globule, having the densest packing density 
of the hexagonal or cubic close packing. When this happens, we discard the individual particles and 
consider instead the globular cluster which they formed. The cluster is porous, so its new radius is 
r = (r/3V2)(3 Ov; /40)18 = (2 /3-V2)(8nv/40)1/3 = 22/3 (nv) 1/3 /(32/827/8). 
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But we do not know the rules by which these particles stick themselves together, whether they 
form a closest-packed globule or some other shapes. It is quite certain that whatever shape they 
are after, they may not retain it for long because there is a limited space within each pore that will 
put constraints on the way they grow. We shall call the growth of clusters into globules mentioned 
above globular formation. 

Other clustering mechanisms possibly include what we shall call the tetrahedra formation. By 
this I mean that each one of our spherical particles attaches itself to three other particles, forming a 
tetrahedron whose side lengths are two times their radius. The next free particle may fit into any of 
the available attachment sites, i.e. the free triangular faces of an existing cluster. We may suppose 
that it always choose the closest one among such sites if this is available; if not, then the next closest 
one and so on. 

These are only two among all the possibilities, namely the globular and tetrahedra formations. 
There could well be others, as well as a mixture of them. On the other hand, each material of 
which the particles are made may decide the particular cluster shape it prefers. Detailed analyses in 
thermodynamics and quantum mechanics are needed if we were to understand this cluster formation 
in continua under spatial constraints. Neither of these is within the scope and time constraint of 
the present work, though both of them merit a detailed investigation which I plan to carry out in 
the future. 

Another problem arises when we come to consider percolation of our membrane. We may, for 
instance, say that it percolates when its structure in three dimensions percolates, or we may say 
that it percolates if there exists a cross section perpendicular to the flow which percolates in two 
dimensions. Since percolation of a cross section implies percolation of the structure but not vice 
versa, these two definitions of percolation due to suspension in membranes are not the same. 

Choosing the percolation of sections as a criterion implies that we consider the superficial 
velocity of the flow whereas choosing the percolation in three dimensions as the criterion means that 
we focus on its interstitial velocity instead. 

The algorithms for the study of percolation by tiny particles due to attrition in membranes 
which proposed here are Algorithm’s 6.19 and 6.20. Algorithm 6.19 prepares the structure while 6.20 
does the percolation simulation. Here both VT and V means the Voronoi tessellation. The appeal 
factor is the probability that a particle will choose to leave a cell via a certain bond. It is weight by 
the gradient of each bond, and is calculated over all bonds going in the downward direction from 
the cell. Transfer grids are square grids which help map the continuous plane at the top layer to 
bonds connected to it, that is to say, it maps a continuous Euclidean plane into discrete grids and 
from there on to bonds. In other words, E? — D? > {b}. 


Algorithm 6.19 Percolation by tiny particles due to attrition in membranes. 


generate a Voronoi tessellation in three dimensions; 
transform the VT into a centroid VT; 
find the cross section of its top layer; 
C+ C(y); 
find transfer grids of C; 
for every cell in VT do 
find the maximum chamber capacity of its cell; 
find appeal factors for all its bonds; 
endfor oD 


Let the gradient of each bond be represented by an angle a that it makes with the horizontal 
plane. Then the gradient can be calculated from the coordinates of the two end points of each bond, 
providing that z2 > 21, from a = tan~!(z12/((Aa)? + (Ay)?)!/2), where the slope is downwards from 
p2 to pi, and as usual zj2 = zg — 21. Algorithm 6.20 describes the membrane percolation simulation 
proposed. 

Algorithm 6.20 Percolation by tiny particles due to attrition in membranes, percolation sim- 
ulation. 


for each time step do 
for all arriving particles do 
find their random arrival position; 
round these positions to the precision of the grids; 
map positions on to bond numbers, using the grids; 
endfor 
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for all particles do 
update distance travelled; 
update chamber crowding; 
find percolation of blocked chambers; 
if chambers percolate then 
terminate the simulation; 
endif 
endfor 
endfor D 


Here we concentrate on the interstitial flow velocity, therefore the percolation is supposed to 
occur when the chambers percolate in three dimensions. Coincidentally, this also makes the calcula- 
tion easier. If we were to choose the percolation of sections as the deciding factor for percolation of 
the membrane, for instance, we would have needed to consider approximately 2n sections in total, 
where n is the number of chambers. With an equal probability for success for all the homogeneous 
sections, this would still leave us on average n sections to consider before we know that a sample 
percolates, if it does, though we would still need to test all the 2n sections in cases where it does 
not. This number 2n arises from the fact that to completely cover all combinations of grouping cells 
into sections we need to consider for each cell two sections for each existing cell, one touching its 
top while the other touches its bottom. 

Sphere packing is a rich field of its own, within which the packing density is generally referred 
to as 7, an efficiency, instead of the usual density symbol p. The rigid packings of spheres vary in 
density from the lowest in loose packing where 7 = 0.06 to the highest, which is shared by the cubic 
and the hexagonal closest packings, 7 = 0.74. A rigid packing is a packing in which all spheres 
touch at least four others, and the points by which each sphere touches its neighbours can neither 
be all in the same hemisphere nor all on an equator, i.e. a greatest circular section. So we can now 
limit the value of p that we shall use to be in accord with 0.06 < 7 < 0.74. In this early stage we 
shall not use a Monte Carlo study to find the probable p, i.e. 7, but will approximate it to be some 
value within the range mentioned. Since the most familiar sphere packing in human history must be 
that by which oranges are stacked at markets, especially open markets like the one at Bolton, which 
gives the efficiency of packing 7 ~ 0.74, we shall assume that this is the way the clusters arrange 
themselves. 

Notice also that piling oranges in a neat tetrahedral shape on a table and a packing them 
into a rectangular box both produce the same crystal structure, that is the face-centred lattice, the 
difference being only in their habits. 

For all intents and purposes the percolation probability of spheres under the influence of the 
van der Waals force must be the same as p, of a face-centred lattice. This is because the biggest 
cluster of both cases will have the same structure and their orientation will determine the orientation 
of the structure. We can do away with the orientation of other minor clusters precisely because they 
are much smaller, which justifies our grossing over their individual shapes and only concern ourselves 
about their statistics, that is to say, their number. I think that stacking oranges into a box is cubic 
close packing while a pile of oranges is hexagonal close packing, but this needs to be checked. 

To find p, of the close-packed cluster of spheres, one needs a program similar to the one 
mentioned in § 6.10, but which would do the job for three dimensions instead of two. For this 
purpose, the program for 2-d tilings mentioned in § 6.10 has been developed further to deal with 
regular lattices in three dimensions. At first I thought that there should be some other way to do 
this instead of having to develop another program for a general lattice in three dimensions, since this 
is already the last week of the project and time is running out. But in the end I found it better to 
spend some time to systematically develop a program for general cases than to opt for some adhoc 
approaches. As a result, a program that creates regular lattices in three dimensions for the purpose 
of percolation study has been written and is listed in § A.30. 

As the 2-d program in § A.6 does for all 2-d regular lattices, this new program can deal with 
all possible lattices in three dimensions. The difficulty is, however, in the meticulous nature of 
identifying all the vertices and links in each unit cell. In this respect, the cubic close packing is 
much simpler to do than the hexagonal close packing. Therefore we shall only do the first one while 
leaving out the second, which ideally could be used for the purpose of comparison. 

Figure 6.11 shows the lattice generated by the program and one which is used for finding the 
percolation threshold. 
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Figure 6.11 Percolation of the cubic close-packed lattice; (a) a unit cell, (b) eight unit cells, 
one from each of the eight groups, (c) network of size5 x 5 x 5 unit cells, which is used in a 
simulation, and (d) the same network with only boundary edges drawn to make it easier to look 
at. 


At present the program only finds py, pe, Ly and x, not Pc, py, Z_e and xj. There may be 
altogether three types of cell and bond pairs, in comparison with two in the 2-d case, depending 
on whether the number of shared vertices required be 3, 2, or 1. For our purpose in the study of 
filtration, we only need to know p,, which, when generated from a 5 x 5 x 5 network as shown in 
Figure 6.11 (c), turns out to be 0.25. All the results from simulations are p, = 0.2501 40.0400, x, = 
8.0645, pe = 0.1320 + 0.0209, x, = 17.1709, while n, = 341 and n, = 1,375. 

This means that when the space will be blocked, i.e. percolates, when it is filled up to one 
quarter of its volume by suspended particles in the form of clusters of the highest packing density. 
Because the cubic close-packed spheres fill 0.74 of the space, this ratio translates into the real volume 
ratio of 0.74 x 0.25 = 0.185, that is 18.5 per cent by volume. If the fluid in our system is water, 
then p = 1,000 kg- m7? and the percentage by volume above is equivalent to a density of of the 
suspended particles of 555 kg -m~?. Notice also in our simulation that the cluster shapes need not 
be convex. 

In the light of the symmetry between particles and space, in other words between particles and 
anti-particles, which give rise to a symmetry and the three types operational regions that I originally 
proposed in a study of traffic congestion (cf § 7), the operational space of our filter may fall into 
three distinct regions when it is subjected to very small particles suspended in a fluid. 

If we specify by p, the ratio of the volume occupied by all the clusters to the total volume, and 
p the density in weight per volume, and if the three regions of operation are labelled I, II and III, 
then in the case of I, 0 < py < 0.25, while for IT, 0.25 < p, < 0.75 and for III, 0.75 < p, <1. In other 
words, for I, II and III, we have respectively 0 < p < 555, 555 < p < 1,665 and 1, 665 < p < 2, 220, 
where the unit of p is kilogram per cubic metre. Generalising this, the regions I, II and III correspond 
respectively to 0 < p< po, = pi, pi < p< Peo = p2 and po < p< pc, where p-, is determined by a 
physical constraint, namely the packing efficiency mentioned. 

Qualitatively speaking, these three regions may correspond to the operational-, blocked and 
non-operational regions. If our system also contains other particles which are larger than the pores, 
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then in Region I filters will operate normally until the effects of blocking or blinding by the large 
particles become prominent, as has been studied in literature (cf Jackson, 1994; Jafferali, 1995). In 
this case fouling of the filter is caused exclusively by the blocking or blinding of these larger particles, 
which result in the formation of cake, and unless p = 0 there will be some blockages of internal pores 
due to the blockage caused by small suspended particles forming cluster. This latter type of blockage, 
which is of our concern, is to be expected due to various reasons. Cluster formation may be caused 
by nonhomogeneity in the concentration of the suspension which raises the concentration in some 
region such that it exceeds p1. 


Additional reduction in the flux is to be expected from two reasons. Firstly blocking clusters 
may block some of the pores. And secondly, the suspended particles together with free clusters, i.e. 
those which are smaller than they could block pores, displace the volume of the liquid surrounding 
them and thereby reduce the flux. The first one of these will produce an effect similar to blinding 
described in literature (cf Jackson, 1994; Jafferali, 1995, ibid.), where no backflushing may recover 
the filters to their virginal state. On the other hand, the majority of those particles and clusters in 
the second scenario is expected to be easily removed when backflushed. Among these latter there 
could yet be some which adhere themselves to the walls, whose fixation defies backflushing. But 
these last ones are expected to be small in number, and thus can be neglected, because we shall 
assume that the combination of p and channelling results in the possibility that pores are blocked 
being very close to either zero or one. 


The channelling effect, or the occurrence of rivulets by some other authors, is the formation of 
preferred paths through a porous media through which liquid and the solids it contains pass. It is 
not yet clear what these rivulets would do to our system. If all the channels channel equally both the 
liquid and the solids, then the blockage along these path may be expected to rise above the average 
value of the whole structure if only because this becomes more probable statistically. But if there 
exist some channels which prefer channelling liquid to particles or vice versa, then the effect they 
produce will vary and become complicated. For example, channels which like to channel particles 
are more likely to find themselves blocked in the end by those particles which pass through them. On 
the other hand those channels which channel liquid better than solids will be less prone to blocking 
on average, but will leave other pores around them with the excess particles, and these latter will 
necessarily become blocked more often than usual. But, for our purpose here, we shall assume that 
it is solid particles that are being channelled. This should raise the possibility of blocking in some 
of the pores by certain amount. Channelling in general needs further investigation which will not 
be covered here. 


Before going on to the next step of our study we should briefly look at the core idea that makes 
the program in § A.30. There are eight types of unit blocks here, compared with the four types in 
the case of the program in § 6.10. These correspond to the area drawn and labelled in Figure 6.12. 
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Figure 6.12 shows the eight areas defined by 


vir VI the eight types of unit blocks they contain. 
/ Area’s from I to IV correspond to those previ- 
7 ously defined for the 2-d program. Although 


unit blocks in the various areas works differ- 
ently, that is to say, they adopts different set 
of vertices from different sources, and create 
different edges, all of them follow the same 
four rules. These four rules in the mnemonic 
not cryptic forms which I use are, take vertices 
l< ; from behind, make front vertices, draw edges 
behind and draw no front edges. With these 
; rules in mind, both programs should become 
| self-illuminating to such extent that no fur- 
\ ther explanation is needed. This set of rules 
II does two things, namely organise vertices and 
then link them with bonds. The unit vec- 
tors in the three directions being orthogonal 
means that we will have a nice and square end 

0 product suitable for a percolation study. 


Figure 6.12 The eight areas defined by the 
eight types of unit blocks. 


> VI 


Pad 


In fact it is wrong to say that all four rules work differently for each unit group. Only the two 
on vertices, viz. the first two, are different. The rest, viz. the last two which concern edges, are the 
same for all basic units. These are the only two crucial tasks with the discovery of which the writing 
of both program becomes worthwhile. 

The input data needs only contain details regarding units in Area II, III and V. Area I, being 
at the origin, is trivial, or should one rather say unique. Area IV can be derived from Area’s II and 
III. Like wise VI is derived from II and V, and VII from III and V. Finally Area VIII turns out to 
be nothing but II, III and V combined. 

So far we have only mentioned the situation where 0 < p< p,. In the second case, where pi < 
p < p2, many more pores are blocked from small particles than in the first case. The concentration is 
already beyond the first critical point. But while the second critical point is still not reached, there 
would still be an infinite cluster of space — in this case the liquid — surrounding the particles. In other 
words the space still percolates. The presence of this infinite cluster, or continuum of the medium, 
means that the filter can still be in operation until it should be blocked or blinded by larger particles 
which individually can physically block the pores as found in existing literature earlier mentioned. 

The third and last case, where po < p < p<, represents the extreme which can be easily 
comprehended. Here the solid particles occupy more space than the liquid does, as a result of which 
the combination is no longer a suspension but a slurry. Clay material produced by this slurry would 
block most of the pores within the structure and make filtration impossible. Backflushing will not 
be effective on filters which have undergone such fate. 


§ 6.11 The second part, flow through the cells 


Next we will investigate briefly the effect that channelling may have on the value of p.. There 
are three cases considered here. The program for this purpose is listed in § A.31. The program 
creates a centroidal Voronoi tessellation in three dimensions. The first set of simulations works on a 
normal case of cell percolation, similar to that in § 6.11 but here the system is a centroidal Voronoi, 
which implies a constraint on cell sizes and distribution. The function perd in it is obtained by 
inputting the Blocked variable instead of generating it internally. 

This section is written along with my developing the programs, so the contents should be easier 
to follow than in other sections. This is not to say that in those other sections I had not kept records 
of things discovered. Every Ph.D. student starts off doing his project knowing that he should write 
along as he goes, and plans to nothing but that. But the truth is that even though we always write, 
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but the way we write develops with our experience. Also, with the the increase in the understanding 
of our problem, we no doubt would be able to write a better description of what we do and how we 
do it. This is unavoidable, and it is probably the reason why we should keep on working. 

So much for an aside. From one hundred generators originally, Voronoi operator is applied 
twice. After the rims has been trimmed there are 280 centroidal Voronoi cells remaining, and this 
is the value of n,. For x, the value is 10.7929, while p, from 2 x 5 simulations is 0.2314 + 0.0602. 

Next investigate the effect of channelling by assuming that the steepest gradient of all the 
bonds arriving at a cell decides how quick it percolates. Here cells are sorted according to their 
steepest gradient of incoming bonds. Working on the same network as previously, if the percolating 
order is such that the steeper the quicker, then p, = 0.2107. But if on the other hand steeper 
incoming bond means slower percolation, then p, = 0.1429. The critical probability is constant in 
this case since the order of percolation is predetermined by the orientation of bonds with respect to 
cells. 

If instead of looking at only a single bond we take the signed summation of bonds entering 
and leaving a cell, then p, = 0.1321 when the criterion is min(}> 6; — >> 6,), and p. = 0.1393 when 
it is max(}> b; — >> b,), where b; and b, are respectively the incoming and outgoing bonds. 

Our studies up to now tell us that if the suspension is homogeneous and there is a rivuleting 
effect, then the location of the blockages made by clusters of suspended particles among all pores of 
the structure is predetermined. This is in contrast with the blinding and blocking of large particles, 
where such location is random. In fact, even for these latter large particles, the location can only be 
random when the particles have a variety of sizes. It can never be wholly random, however, if this 
is not the case, since in the former case the randomness is introduced by the distribution of sizes 
which is random, but in the latter the blinding or blocking will be determined by the size of pore 
openings which is fixed by the geometry of each network. The randomness then can only be in the 
order not location of blockings. 

The next study is a combination between continuum and network percolations. Here suspended 
particles are grouped into quanta, each channelling through a path or rivulet of interstitial distance 
with some interstitial velocity. 

At the top, the layer of the structure where the incoming particles arrive, is a cross section 
all the cells of which are gridded to provide means to determine the path at the beginning of each 
quantum. The partitions in this layer is conveniently found by cutting some faces of the convex hull 
of each cell in the top layer by the plane z = 0.9zmax, where Zmax is the maximum z-coordinate of 
all the cells under consideration. 

Given coordinates of two points, (#1, 41, 21) and (22, y2, 22), and a plane equation z — a = 0, 
we may think of the plane equation as being one coordinate given, z = a, and find the coordinates 
of intersection between a line passing through the two points and the plane from the parametric 
equations for the line. Parametric equations are in fact interpolation done on each of the coordinates. 
In this case, which is useful when finding the intersection between an edge of a triangle and a plane 
perpendicular to some coordinate axis, the parametric equations are x = 21 + Zot, y = yi + Yrat 
and z = 21 + 212t. From the plane equation z = a, therefore t = (a — 21) /z12 = (a — 21)/(z2 — 21). 


The program being written finds the inter- i 


section of cells with the horizontal plane by —to9 
finding the intersection of the faces of its con- 
vex hull with the same. Since every face of the 
convex hull is a triangle, the program essen- 
tially finds intersection between edges of these 
triangles and the horizontal plane. The result 400; 
obtained from an intermediate state during 
the course of development of the program is 
shown in Figure 6.13. The partitions look in- 600/ 
complete because the picture is taken as a test 


while developing the program as mentioned. 7 
A picture with the same degree ofincomplete- gu) 
ness as this one is not to be obtainable from 

the completed program. 900; 


Figure 6.13 Intersection of convex hull faces 999 
0 


and the horizontal plane. 200 400 600 800 1000 


nz = 550 
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For each face of the convex hull that intersects the plane, there will be two points of intersection 
arising from the two edges of the triangle intersecting it. Let (v1, y1) and (x2, y2) represent these 
two points. Then we may scan up in the y direction finding x = a + (y — y1)@12/y12 for each y 
along the way, and then scan in the x direction, this time finding instead y = y1 + ( — 21) y12/212. 
Afterwards we could fill in the space by scanning along the x direction for all y positions. 

Figure 6.14 shows the progress of my programming the codes, step by step, trying to close 
all the partitions such that no gaps remain. Two problems have been discovered, namely those of 
rounding and precision. Before the correction the resul looks like Figure 6.14 (a) and (b), and after 
rounding problem corrected Figure 6.14 (c). 
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(a) (b) 
Figure 6.14 Correcting the effect of rounding, in other word discretisation. (a) The gaps 
resulted from rounding or discretisation, (b) a closed-up view of (a) and (c) partial remedy 
where roundings have been solved but with the degree of precision not yet raised. 


After having corrected the problem regarding precision, by increasing the number of steps 
when calculating x or y, the result still misses several walls, the cause of which is still unknown at 
present. This is shown in Figure 6.15. 


0 


ase 0 50 100 150 200 

nz = 1051 

(b) 
Figure 6.15 After having increased the precision to ten times the previous value, (a) and a 
closed-up view, (b). 


It seemed at first that there might be some triangular faces missing from the surface of the 
convex hull, which would have accounted for the missing boundary in the section. But after having 
tested the minimum number that an edge of each hull appears as the edges of all its triangular faces, 
and see that it’s value is correctly two, this becomes out of question. 

The figures, viz. Figure 6.13, 6.15 and 6.15, are produced from the spy command in Matlab, as 
a result of which the z-axis runs downwards while the y-axis runs to the right. This command looks 
at a matrix from above as we look at a map. In the present case our matrix is a full-, not sparse 
matrix. The number written at the bottom is the number of all its nonzero components, which is 
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less than the number of times that we calculate them since we need to calculate some of the points 
more than once in order to increase the precision to eliminate gaps in other places. 


The command spy is used more often with sparse matrices since these are often too large to 
list, and listing their members in pairs makes it difficult to visualise. As an example, Figure 6.16 
(a) is what we get when we spy our neighbour matrix necc, while in Figure 6.16 (b) are all the 
neighbours that the 100 cell has. Cells which have few neighbours generally live along the border. 


For example the 110* cell has only three neighbours, and it is located not far from the lower zx limit. 
OF nr Mee be Tee 
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(b) 


Figure 6.16 The result when we spy our neighbour matriz. Here (a) shows the neighbours of 
every cells while (b) only shows those of the 100 cell. 


At this stage the neighbour matrix of the program is double-checked, and find that it includes 
poorly defined neighbours, that is those which only have one vertex in common. Therefore the 
program is first altered to make it only look for neighbours who share at least two vertices. But this 
has shown no noticeable changes in the results that we have so far. 


The next step is to colour the cells. Like a child painting a picture, there are so many ways 
one can paint or label the tiles of a tessellation. For example one could draw a vertical line first 
and then branch out either horizontally or diagonally. On the other hand, one could also draw the 
diagonals first Yet another way is to expand radially, spiraling outwards. With parallel computing 
we could also divide the area into domains, paint each domain, and then merge the resulted areas 
together. 


But here we opt for painting the cells by scanning horizontally, moving upwards in layers. 
Once past a wall, the program moves on until the next wall is reached while gathering all the grids 
that are between the two walls into one group. It then colour the whole group by a colour picked 
up from one layer below it. If no colours exist, then it creates a new colour which in turn gradually 
propagates upwards this way until the upper wall is reached. 
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The intersection of bonds will not always coincide 
with the intersection of the cell, neither is the pro- 
jection of a cell perpendicular to the plane the sec- 
tional plane of the cell. This is because of the three 
possible situations shown in Figure 6.17. In Fig- 
ure 6.17 (a) the cell section contains several points, 
while both (b) and (c) contain none. 


(c) 


Figure 6.17 Sections of bonds and cells. 


And here sadly the time runs out, so I will suffice myself to describing what I see should be 
done next. Up to now we have a three-dimensional network and its top section. We also have the 
list of all its bonds, which contains the connections and the draining angles sorted in a descending 
order. 

Next we should find a mapping from each cell section to the corresponding nuclei. Then the 
times it takes to traverse each bond must be calculated. This time for each bond is then divided by 
half, one belonging to each of the two cells connected by the bond. 

When it comes to bombarding the filter with our small particles, we can not keep track of 
millions of particles and therefore we should quantise them into units. These units or quanta can 
then be treated as individual particles. When a quantum enters a cell, it is assigned ¢, the time to 
reach the nucleus. Later time sees this ¢ decreases in steps until it finally reaches its destination, the 
nucleus. Once there, it is assigned the next bond to go along, taking into account what bonds are 
available at the time and their comparative probabilities, which in turn depend on their gradient as 
mentioned. When this is decided, it is given ¢;, the time it would take to reach the border that lies 
at mid point of the bond. 

This goes on forever, apart from that at each time step we look to see whether the blockages 
in our filter has percolated. After updating the list of blocked cells, if we find that percolation has 
occurred then the simulation would end. Percolation occurs in each cell whenever its concentration 
has reached a certain value. This value we have found in § 6.11 to be 18.5 per cent by volume. 

To calculate the flux decrease we find instead the decrease in the superficial area. This is 
calculated from the total volume of the void subtracted by the volume of all cells that had percolated, 
and then subtracted by the total volume of solid particles which are suspended inside the network. 
The area of the cross section is then the volume which remains divided by the thickness of the filter. 

It is not a little to have to leave things unfinished after having started it off. But as one New 
Zealander poet says, ‘Alone we are born and die alone, yet see the red-gold cirrus over the snow 
mountains shines. Upon the up-land road ride easy, stranger. Surrender to the sky your heart of 
anger.’ And with this we go on to the next chapter. 
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§ 7. Percolation in traffic modelling 

It does not take much imagination for an average person to see that traffic congestion is a 
percolative process. As I have lived and studied in Bangkok where, at that time, the traffic jams 
were renowned. Now I relate traffic jam to percolation of the network of streets. I remember once 
in 1990 it took me more than five hours on a bus to travel the distance of ten kilometers along the 
Sukhumvit Road in Bangkok. The newspaper of the following day reported that that traffic jam 
had been on the BBC news broadcasted all over the world. 

Soon after having decided to investigate the application of percolation to traffic, I gave a 
poster presentation at the Fourth Annual Conference of Thai Researchers in Japan in 1997, the 
proceedings of the event of which also contain a number of abstracts from some of the researches 
that I was currently working on (Tiyapan, 1997, KNT5(i) to KNTs(iv)). The poster presentation was 
considerably a success. 

The paper listed in § E.20 of Tiyapan (2003, KNTs(ii)) is a reproduction of what I had submitted 
to the Journal of Statistical Physics in 1997. I tried but could not receive any reply for the paper I 
submitted. Afraid that the ideas in it could have been plagiarised, I decided to put it in the book 
above quoted. 

In 2002 I did some more works on it. I have developed the theory further and done a number 
of simulations. I plan to carry on doing research along this line in the future, together with another 
area where percolation is used to explained economics transition (Tiyapan, 1997, KNTs(viii)). 

In his paper submitted to the Journal of Statistical Physics Tiyapan (1997) introduces a new 
idea of considering the development of clusters in both phases at the same time. Applied to the 
context of traffic network, these phases are namely cars and spaces. Furthermore, because both 
phases reside in one and the same network, there is a symmetry which divides the probability 
space into three regions, symmetric with respect to p = 0.5. This helps divide the traffic condition 
into three regions as existing literature in traffic study at the time described, namely free flowing, 
congested and stand still. In particular, this idea explains the difference between the congested and 
the stand still states. There has been no reply from the journal. 


Tn the 2002 simulation the networks are drawn 

which have their vertices as points where two _ 

or more roads meet one another. The pro- | 

posed study is to compare the robustness of space ears 


two traffical networks by comparing their per- 
colation thresholds. When a new motorway 
is planned, for example a ring road around a 
city, the two networks, one with the ring road 
and the other one without, can be simulated 
to find their percolation thresholds and then 
these values compared. 

Traffic status or condition when the per- 
colation probability p, of cars is more than 


I 
0 (1-Pc) 0.5 Pc 1 


0.5. The traffic condition is generally described 

as free-flowing, congested, or stand-still. For 

cases where p, > 0.5 as the one shown in Fig- free-flowing congested stand-still 
ure 7.1 the stand-still traffic corresponds to x >\« _ >| 


the situation where only cars have percolated 
but not the space, that is the places on the 
road available and accessible to the cars. 


Figure 7.1 Traffic status, p. > 0.5. 


A free-flowing traffic is where only the space but not the cars has percolated, and a congested 
traffic is that when neither of the two has. Traffic status or condition when p, is less than 0.5. 
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On the other hand, in the case where p, < 0.5 
shown in Figure 7.2 the critical probability 
of cars is now on the left hand side of the 
middle line instead of on the right hand side 
thereof in the previous case where p, > 0.5. 
Definitions of the free-flowing and the stand- 
still statuses remain the same, but the con- 
gested traffic is now the traffic where both 
cars and space have percolated. An interest- 
ing question is whether or not there is a differ- 
ence between the congested area in Figure 7.1 
0 Pc 05 (1-Pe) 1 and the one in Figure 7.2. A modelling algo- 
rithm which studies the percolation of a traffi- 

cal networks starts by finding all the vertices 

free—flowing congested stand-still and edges forming a network from the road 
|~« ~~ oe = data. Then a blocking algorithm operates by 
randomly shutting off one edge after another 

until the network percolates when the critical 
probability of the network may be calculated. 


Figure 7.2 Traffic status, p. < 0.5. 


A control algorithm for the real time traffic control, however, is as shown in Figure 7.4. 


Let C means cars have percolated and S 
means space has percolated, then in Figure 
7.4 the cases I, IT, and III are respectively 
AC AS, (AC A7S)V(CAS), and (CA-S). 
Case I is the normal congestion, nothing 
to worry about. Examples of the control 
schemes used in Case II are overriding of the 
find percolations traffic lights by the manual control of traffic 
of cars and space polices at certain strategic points, tempo- 
rary one-way systems, bird-eye view obser- 
vation and feedback from helicopters, traf- 


case fic control centre, distributed control centres 

i together with traffic radio channel broad- 

Y cast. And lastly the most important and 

cBiitio laches eaiescuey critical Case III the emergency plan of which 

igi aaa applied measure may include directing all cars away from 


congested clusters or if necessary out from 
the city, and directing all incoming traffic 
such that no more cars may enter the city 
until the emergency status ends. 


Figure 7.4 Proposed traffic control in real 
time.. 


A congested cluster can be broken up by forming a one way flow channel cutting through it 
which leads cars away from the cluster. How far the channel needs to go before letting the cars 
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on it circle and seep back into town depends on how badly congested the traffic is. Although it 
is a normal practice to lead cars along a long detour because this brings in more road surface and 
thus enlarges the network, identifying the percolating cluster and cutting it into two or more parts 
by guarded flow channels is much less common or even unheard of. The latter seems to be more 
important and will lead to a better and more effective control, namely the control of the percolating, 
in other words the biggest, cluster. 

The percolation probability is important for the networks of traffic both inside cities and among 
them. It shows the degree of connectivity of the area being considered. Urban road networks have a 
general character which differ from one country to another, the simplest construction of which seems 
to be that of the square lattice. One can find this theme of the square lattice and its variations, 
similar to the dislocations and defects found in minerals, in America. Examples are Denver, Aspen, 
Durango, Pueblo, Salida and La Junta in Colorado; Boise, Pocatello and Twin Falls in Idaho; Butte, 
Bozeman, Coeur d’Alene, Kalispell in Montana; and Cheyenne, Laramie and Sheridan in Wyoming 
(cf Florence et al). The Great Junction in Colorado and the Great Falls in Montana are very close 
to being perfect square lattices. As more examples of these (cf Collins USA, 1999), in Arizona there 
are Phoenix, Yuma, Tucson; in California Bakersfield, Central San Francisco, Central Sandiego, 
Fresno, (Central) Los Angeles and vicinity, Modesto, Sacramento; in Colorado Fort Collins, Denver 
and vicinity, Greeley; in Florida Central Miami; in Georgia Central Atlanta; in Illinois Champaign 
and Urbana, (Central) Chicago and vicinity,, Quad Cities, Rockford; in Indiana Fort Wayne and 
Indianapolis; in Kansas Topeka and Wichita; in Louisiana Central New Orleans; in Maryland Central 
Baltimore; in Minnesota Central Minneapolis and Central St. Paul; in Missouri Central Kansas City; 
in Nebraska Lincoln; in Nevada Las Vegas; in New York Manhattan; in Oklahoma Lawton, Norman, 
Oklahoma City and Tulsa; in Pennsylvania Central Philadelphia; in South Dakota Sioux Falls; in 
Amarillo, Central Houston and Lubbock; in Utah Central Salt Lake City; in Washington Central 
Washington D. C. and Central Seattle; and in Wisconsin there is Central Milwaukee. One example 
in Canada is Toronto. The square lattices of these cities are sometimes cut through by motorways 
or interstate highways as is the case in Amarillo, Texas. Or they can be surrounded by a ring road 
or a county highway as is what happens with Lubock, also in Texas. 


When the percolation probability is greater than 0.5, we have the interval p. + (pe — 0.5) where 
neither the blocked nor the free roads percolate. If this interval is narrow, that is if p. — 0.5 is small, 
then within this interval the condition of the traffic is very sensitive, and even a seemingly small 
change may lead to a standstill or instead to a free-flowing traffic. This is easily visualised, since in 
such situation there would be small islands of free-flowing roads within a large congested cluster, and 
vice versa small clusters of congested roads within an otherwise noncongested area. Despite their 
sizes, such small islands of anomaly in either of the phases are particularly important. Moreover, 
their importance increases the closer p, is to 0.5. The same characteristic happens in politics where, 
in the case of coalition governments, a minority party which has relatively few representatives can 
become critically important and influential to the major party when the latter needs them in order 
to be able to govern (Ireland, 2002). 

Plan to build a ring road around a city usually includes flyovers, overbridges, and tunnels 
in order to avoid having intersections. The design philosophy used in Europe is to have heavier 
traffic goes under a lighter one, which results in either the ring road going into a tunnel or having 
overbridges or viaducts over it. The philosophy used in Thailand which used to be for the heavier 
traffic go over the lighter traffic, which puts a more severe limit on the weight of trucks and is 
therefore not economic in the long run, but this has started to changed, if only to follow the practice 
of the west. Ring roads do not necessarily resemble a circle, as in the case of the circular speedway 
proposed around Saint Petersburg (Petersburg, 2001), which is in the shape of a cashew nut. It will 
link the arterial roads of the city to Helsinki, Kiev, Moscow, Murmansk, and Tallinn into a network. 
The route proposed is 155 kilometres in length, has 31 bridges, 16 overbridges, 55 viaducts, and 
will support the volume of traffic of 21 million tonnes. The implementation and contract work is 
looked after by a joint-stock company KAD Sankt-Peterburga, under the order of St. Petersburg 
and Leningradskaya Oblast. 

One classical example of a ring road is the M25 motorway which forms a circumscribed ring 
around London. Another more recent example is the M60 orbital motorway around Manchester. In 
the case of M60, various sections of existing motorways have been put together and renumbered. 
The northwest quarter used to be M62, the southwest one M63, and parts of the remaining used 
to be M66. The motorway forms a complete ring around Manchester since 2000 with the opening 
of the final northeastern part which stretches from Denton to Prestwich. There is another smaller 
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ring, an inner ring, which is formed by the A6010 and A576. Orbital motorways around cities have 
now become indispensable and are the hallmark of a city. 


§ 7.1 Percolation of road networks 


I began my study of the traffic networks in 1997 while in Japan. In the same year I submitted 
two papers, 

Tiyapan in 1997 (KNTs5(ix)) has a novel idea of linking percolation to the management of street 
congestion. In his work submitted that year to the Journal of Statistical Physics, also in § E.19 
and § E.20 of Tiyapan (2003, KNTs(iii)), electronically though the paper became lost and was never 
published, he explains how three phases could result from structures containing only two phases. 
Recently there have been some papers published in 2001 in the Journal of Statistical Physics which 
used this very idea that Tiyapan introduced in 1997. But what has passed is past and past cure, we 
need to think about it no further now. Everything comes from and goes to God, let that suffice! I 
have developed a program and a procedure for finding the percolation thresholds of road networks. 
These programs are listed in § A.27. They are different from the usual percolation programs used for 
other kinds of networks including Voronoi. Because of the existence of flyovers and elevated express 
ways, we can not use the duality operator to transform a network of vertices and edges into one 
of cells and bonds. Having said that, the said transformation could become useful in the future in 
some other applications, in some other areas or even within the study of traffic network itself. But 
at present I have only one application in mind for the dual networks of roads, and that is related to 
fire prevention where such an application is by no mean obvious. Therefore, here the boundary is 
defined again to obtain the cells which represents the zones. 

The codes also contain data of several towns and cities, namely Amsterdam, Brussel, Freiburg 
and Manchester. There are three main datasets. The first one contains a list of the coordinates of 
all the vertices. The second one is a list of edges, together with the numbers of the cells that each 
of them connects. The third one contains the coordinates of turns in each of the roads listed in the 
second dataset. This is in order for the graph to look like the actual roads it represents, instead of 
containing only straight lines, as would have been the case were the windings of the roads not to be 
taken into account. Also, these coordinates will make it possible to calculate the true length of each 
road. Even though we have no use for these lengths at the present stage, future developments may 
need them. 

The procedures developed for gathering the data and processing them can be carried out by 
a single researcher, and require no sophisticated tools. Were these tools become available in the 
future, the former could be adjusted to accommodate them. 

To my surprise, so much so that I at first thought that there was something wrong with the 
program, the mean coordination number of the road networks of Manchester turns out to be exactly 
3. The second simulation gives x, = 3.0513, which is still very close to three. 

Mancheter Mancheter (fire control area) 


1000 metres ——————__ 1000 metres 
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(a) 


(b) 


Figure 7.5 For Manchester in this picture, ny = 220, ne = 380, ne = 25, ny = 50, Ly = 
8.0000, %- = 4.7278, t = 4.0000 and x = 6.5200. 


Figure 7.6 shows the largest clusters pl 


otted against p, the probability that each vertex, edge, 


cell or bond respectively for (a), (b), (c) or (d). Here p® is the percolation probability of the 
space in the network of vertices, and similarly for p%, p$ and p?. Each plot represents one of 


the runs of simulation which, for networks 
critical probabilities all over the place. But 


of these sizes and variances, literally distributes the 
the plots of the largest cluster sizes always look very 


symmetrical. This seems to suggest that these sizes may represent the point of percolation better 


than the percolation probability. 
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Figure 7.6 Manchester. Plots of the largest clusters of, (a) vertices, (b) edges, (c) cells and 


(d). Here (py, p§) = (0.5455, 0.7182), (pe 


,p’) = (0.6152, 0.4939), (De,p8) = (0.6400, 0.6000) 


and (pp, p§) = (0.1200,0.5000) respectively for the cases of (a), (b), (c) and (d). 


The mean percolation probabilities from (2x10) similar simulations as the two shown in each of 


the four cases of Figure 7.6 are py = 0.6723 4 
and py = Er[0.3720,0.1149]. 


| 0.0762, p. = Er[0.6083,0.0838], p. = Er[0.6240,0.0921] 


Next consider Amsterdam in The Netherlands. The map for our purpose is shown in Figure 
7.7. Areas shown in Figure 7.7 (b) are usually bound by main roads or the rims of the picture. 


There is no definite relations between vertic 
other. 


es and edges on one hand, and cells and bonds on the 
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Amsterdam Amsterdam (fire control area) 
ae, 
LE 
1000 metres 1000 metres 
(a) (b) 


Figure 7.7 Amsterdam in this picture has n, = 487, ne = 745, ne = 18, np = 35, Ly = 8.0518, 
Le = 4.7007, &. = 3.8889 and x, = 6.8000. 
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Figure 7.8 Amsterdam. Plots of the largest clusters of, (a) vertices, (b) edges, (c) cells and 
(d). Here (p,,p®) = (0.7187,0.7146), (pe,p2) = (0.6309,0.5597), (pe,p3) = (0.7778, 0.4444) 
and (p5,P,) = (0.4286, 0.2571) respectively for the cases of (a), (b), (c) and (d). 


From (2 x 10) simulations we obtain p, = 0.7374 + 0.0500 and p, = 0.6328 + 0.0595, whereas 
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from (2 x 11) simulations p, = 0.5960 + 0.1395 and jp 0.3922 + 0.1395. 


Then consider the road networks of Brussel in Belgium, as shown in Figure 7.9, to gether with 
the fire control area, and the largest cluster sizes in Figure 7.10. Fire control strategy is only one 
possible application to which the percolation of areas. Many other applications which are similar 
in nature, for instance emergency evacuation zones, earthquake evacuation zones, zones prepared as 
measure against a terrorist gas attack, etc. There are also other applications, for example strategic 
areas in market planning and the study of mineral deposits. 

Brussels Brussels (fire control area) 
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Figure 7.9 Brussel has ny = 287, ne = 454, Ne = 21, np = 42, Ly = 3.1568, xe = 8.1568, 
Leo = 4.0000 and xp = 7.0476. 
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(c) (d) 


Figure 7.10 Brussel. The largest clusters when percolate by means of (a) vertices, (b) edges, (c) 
cells and (d). Here (py, p’) = (0.6167, 0.6411), (pe, p?) = (0.6057, 0.6586), (pe, p2) = (0.6667, 
0.5714) and (py, pf) = (0.3095,0.4286) respectively for the cases of (a), (b), (c) and (d). 


From (2 x 10) simulations on Brussel we obtain p, = 0.6580 + 0.0771, pe = 0.6205 + 0.0540, 
Pe = 0.6286 + 0.1043 and pp, = 0.3786 + 0.0924. 


And then consider a small town Freiburg in Germany, where the many roads that are reserved 
for pedestrians only seem at a first glance to have altered much of the structure. But simulations 
have shown that the percolation probabilities remain comparable with networks of other towns. The 
area- and bond coordination numbers obtained for Freiburg are rather low compared with other 
towns. This could mean that the emergency properties of the town is different from those of others. 
Its lower connectivity could mean that it is more robust than others against an attack or in the 
face of catastrophe. But it could also make it more difficult to evacuate from an area. More precise 
relationship between the valence and the interpretation in terms of physical networks can only be 
possible by more extensive investigations in the future. 


Freiburg Freiburg (fire control area) 
500 metres 500 metres 
(a) (b) 


Figure 7.11 Freiburg has ny = 75, Ne = 102, ne = 10, ny = 15, Ly = 2.7200, x = 4.4510, 
te = 8.0000 and xp = 4.5338. 
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